source: codes/icosagcm/branches/SATURN_DYNAMICO/ICOSAGCM/src/caldyn_gcm.f90 @ 314

Last change on this file since 314 was 298, checked in by milmd, 10 years ago

Less output messages are written. On 20000 cores it is better. In LMDZ, only master of MPI and OpenMP can write.
File size: 39.5 KB
Line 
1MODULE caldyn_gcm_mod
2  USE icosa
3  USE transfert_mod
4  PRIVATE
5
6  INTEGER, PARAMETER :: energy=1, enstrophy=2
7  TYPE(t_field),POINTER :: f_out_u(:), f_qu(:), f_qv(:)
8  REAL(rstd),SAVE,POINTER :: out_u(:,:), p(:,:), qu(:,:)
9  !$OMP THREADPRIVATE(out_u, p, qu)
10
11  TYPE(t_field),POINTER :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:)
12  TYPE(t_field),POINTER :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
13
14! temporary shared variable for caldyn
15  TYPE(t_field),POINTER :: f_theta(:)
16  TYPE(t_field),POINTER :: f_pk(:)
17  TYPE(t_field),POINTER :: f_geopot(:)
18  TYPE(t_field),POINTER :: f_wwuu(:)
19  TYPE(t_field),POINTER :: f_planetvel(:)
20   
21  INTEGER :: caldyn_conserv
22 !$OMP THREADPRIVATE(caldyn_conserv)
23 
24  TYPE(t_message) :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu
25
26  PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields, &
27       req_ps, req_mass
28
29CONTAINS
30 
31  SUBROUTINE init_caldyn
32    USE icosa
33    USE exner_mod
34    USE mpipara
35    USE omp_para, ONLY: omp_master
36    IMPLICIT NONE
37    CHARACTER(len=255) :: def
38    INTEGER            :: ind
39    REAL(rstd),POINTER :: planetvel(:)
40 
41    def='energy'
42    CALL getin('caldyn_conserv',def)
43    SELECT CASE(TRIM(def))
44    CASE('energy')
45       caldyn_conserv=energy
46    CASE('enstrophy')
47       caldyn_conserv=enstrophy
48    CASE DEFAULT
49       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', &
50            TRIM(def),'> options are <energy>, <enstrophy>'
51       STOP
52    END SELECT
53    IF (is_mpi_root .AND. omp_master) PRINT *, 'caldyn_conserv=',def
54
55    CALL allocate_caldyn
56
57    DO ind=1,ndomain
58       IF (.NOT. assigned_domain(ind)) CYCLE
59       CALL swap_dimensions(ind)
60       CALL swap_geometry(ind)
61       planetvel=f_planetvel(ind)
62       CALL compute_planetvel(planetvel)
63    END DO
64
65  END SUBROUTINE init_caldyn
66
67  SUBROUTINE allocate_caldyn
68  USE icosa
69  IMPLICIT NONE
70
71    CALL allocate_field(f_out_u,field_u,type_real,llm) 
72    CALL allocate_field(f_qu,field_u,type_real,llm) 
73    CALL allocate_field(f_qv,field_z,type_real,llm) 
74 
75    CALL allocate_field(f_buf_i,   field_t,type_real,llm)
76    CALL allocate_field(f_buf_p,   field_t,type_real,llm+1) 
77    CALL allocate_field(f_buf_u3d, field_t,type_real,3,llm)  ! 3D vel at cell centers
78    CALL allocate_field(f_buf_ulon,field_t,type_real,llm)
79    CALL allocate_field(f_buf_ulat,field_t,type_real,llm)
80    CALL allocate_field(f_buf_v,   field_z,type_real,llm)
81    CALL allocate_field(f_buf_s,   field_t,type_real)
82
83    CALL allocate_field(f_theta, field_t,type_real,llm,  name='theta')   ! potential temperature
84    CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk')
85    CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot')  ! geopotential
86    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu')
87    CALL allocate_field(f_planetvel,  field_u,type_real, name='planetvel') ! planetary velocity at r=a
88
89  END SUBROUTINE allocate_caldyn
90
91  SUBROUTINE caldyn_BC(f_phis, f_wflux)
92    USE icosa
93    USE mpipara
94    USE omp_para
95    IMPLICIT NONE
96    TYPE(t_field),POINTER :: f_phis(:)
97    TYPE(t_field),POINTER :: f_wflux(:)
98    REAL(rstd),POINTER  :: phis(:)
99    REAL(rstd),POINTER  :: wflux(:,:)
100    REAL(rstd),POINTER  :: geopot(:,:)
101    REAL(rstd),POINTER  :: wwuu(:,:)
102
103    INTEGER :: ind,i,j,ij,l
104
105    IF (omp_first) THEN
106       DO ind=1,ndomain
107          IF (.NOT. assigned_domain(ind)) CYCLE
108          CALL swap_dimensions(ind)
109          CALL swap_geometry(ind)
110          geopot=f_geopot(ind)
111          phis=f_phis(ind)
112          wflux=f_wflux(ind)
113          wwuu=f_wwuu(ind)
114         
115          DO ij=ij_begin_ext,ij_end_ext
116              ! lower BCs : geopot=phis, wflux=0, wwuu=0
117              geopot(ij,1) = phis(ij)
118              wflux(ij,1) = 0.
119              wwuu(ij+u_right,1)=0   
120              wwuu(ij+u_lup,1)=0   
121              wwuu(ij+u_ldown,1)=0
122              ! top BCs : wflux=0, wwuu=0
123              wflux(ij,llm+1)  = 0.
124              wwuu(ij+u_right,llm+1)=0   
125              wwuu(ij+u_lup,llm+1)=0   
126              wwuu(ij+u_ldown,llm+1)=0
127          ENDDO
128       END DO
129    ENDIF
130
131!    !$OMP BARRIER
132  END SUBROUTINE caldyn_BC
133   
134  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
135       f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
136    USE icosa
137    USE disvert_mod, ONLY : caldyn_eta, eta_mass
138    USE vorticity_mod
139    USE kinetic_mod
140    USE theta2theta_rhodz_mod
141    USE wind_mod
142    USE mpipara
143    USE trace
144    USE omp_para
145    USE output_field_mod
146    IMPLICIT NONE
147    LOGICAL,INTENT(IN)    :: write_out
148    TYPE(t_field),POINTER :: f_phis(:)
149    TYPE(t_field),POINTER :: f_ps(:)
150    TYPE(t_field),POINTER :: f_mass(:)
151    TYPE(t_field),POINTER :: f_theta_rhodz(:)
152    TYPE(t_field),POINTER :: f_u(:)
153    TYPE(t_field),POINTER :: f_q(:)
154    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
155    TYPE(t_field),POINTER :: f_dps(:)
156    TYPE(t_field),POINTER :: f_dmass(:)
157    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
158    TYPE(t_field),POINTER :: f_du(:)
159   
160    REAL(rstd),POINTER :: ps(:), dps(:)
161    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:), dtheta_rhodz(:,:)
162    REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
163    REAL(rstd),POINTER :: qu(:,:)
164    REAL(rstd),POINTER :: qv(:,:)
165
166! temporary shared variable
167    REAL(rstd),POINTER  :: theta(:,:) 
168    REAL(rstd),POINTER  :: pk(:,:)
169    REAL(rstd),POINTER  :: geopot(:,:)
170    REAL(rstd),POINTER  :: convm(:,:) 
171    REAL(rstd),POINTER  :: wwuu(:,:)
172       
173    INTEGER :: ind
174    LOGICAL,SAVE :: first=.TRUE.
175!$OMP THREADPRIVATE(first)
176   
177    ! MPI messages need to be sent at first call to caldyn
178    ! This is needed only once : the next ones will be sent by timeloop
179    IF (first) THEN
180      first=.FALSE.
181      IF(caldyn_eta==eta_mass) THEN
182         CALL init_message(f_ps,req_i1,req_ps)
183      ELSE
184         CALL init_message(f_mass,req_i1,req_mass)
185      END IF
186      CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
187      CALL init_message(f_u,req_e1_vect,req_u)
188      CALL init_message(f_qu,req_e1_scal,req_qu)
189!      IF(caldyn_eta==eta_mass) THEN
190!         CALL send_message(f_ps,req_ps)
191!         CALL wait_message(req_ps) 
192!      ELSE
193!         CALL send_message(f_mass,req_mass)
194!         CALL wait_message(req_mass) 
195!      END IF
196    ENDIF
197   
198    CALL trace_start("caldyn")
199
200      IF(caldyn_eta==eta_mass) THEN
201         CALL send_message(f_ps,req_ps) 
202         CALL wait_message(req_ps) 
203      ELSE
204         CALL send_message(f_mass,req_mass) 
205         CALL wait_message(req_mass) 
206      END IF
207
208    CALL send_message(f_u,req_u)
209    CALL wait_message(req_u)
210    CALL send_message(f_theta_rhodz,req_theta_rhodz) 
211    CALL wait_message(req_theta_rhodz) 
212   
213!    CALL wait_message(req_u)
214!    CALL wait_message(req_theta_rhodz)
215
216    SELECT CASE(caldyn_conserv)
217    CASE(energy) ! energy-conserving
218       DO ind=1,ndomain
219          IF (.NOT. assigned_domain(ind)) CYCLE
220          CALL swap_dimensions(ind)
221          CALL swap_geometry(ind)
222          ps=f_ps(ind)
223          u=f_u(ind)
224          theta_rhodz = f_theta_rhodz(ind)
225          mass=f_mass(ind)
226          theta = f_theta(ind)
227          qu=f_qu(ind)
228          qv=f_qv(ind)
229          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv)
230       ENDDO
231
232       CALL send_message(f_qu,req_qu)
233       CALL wait_message(req_qu)
234
235       DO ind=1,ndomain
236          IF (.NOT. assigned_domain(ind)) CYCLE
237          CALL swap_dimensions(ind)
238          CALL swap_geometry(ind)
239          ps=f_ps(ind)
240          u=f_u(ind)
241          theta_rhodz=f_theta_rhodz(ind)
242          mass=f_mass(ind)
243          theta = f_theta(ind)
244          qu=f_qu(ind)
245          pk = f_pk(ind)
246          geopot = f_geopot(ind) 
247          CALL compute_geopot(ps,mass,theta, pk,geopot)
248          hflux=f_hflux(ind)
249          convm = f_dmass(ind)
250          dtheta_rhodz=f_dtheta_rhodz(ind)
251          du=f_du(ind)
252          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du)
253          IF(caldyn_eta==eta_mass) THEN
254             wflux=f_wflux(ind)
255             wwuu=f_wwuu(ind)
256             dps=f_dps(ind)
257             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
258          END IF
259       ENDDO       
260       
261    CASE(enstrophy) ! enstrophy-conserving
262       DO ind=1,ndomain
263          IF (.NOT. assigned_domain(ind)) CYCLE
264          CALL swap_dimensions(ind)
265          CALL swap_geometry(ind)
266          ps=f_ps(ind)
267          u=f_u(ind)
268          theta_rhodz=f_theta_rhodz(ind)
269          mass=f_mass(ind)
270          theta = f_theta(ind)
271          qu=f_qu(ind)
272          qv=f_qv(ind)
273          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv)
274          pk = f_pk(ind)
275          geopot = f_geopot(ind) 
276          CALL compute_geopot(ps,mass,theta, pk,geopot)
277          hflux=f_hflux(ind)
278          convm = f_dmass(ind)
279          dtheta_rhodz=f_dtheta_rhodz(ind)
280          du=f_du(ind)
281          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du)
282          IF(caldyn_eta==eta_mass) THEN
283             wflux=f_wflux(ind)
284             wwuu=f_wwuu(ind)
285             dps=f_dps(ind)
286             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
287          END IF
288       ENDDO
289       
290    CASE DEFAULT
291       STOP
292    END SELECT
293
294!!$OMP BARRIER
295    IF (write_out) THEN
296
297!       IF (is_mpi_root) PRINT *,'CALL write_output_fields'
298
299! ---> for openMP test to fix later
300!       CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
301!            f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p)
302       CALL un2ulonlat(f_u, f_buf_ulon, f_buf_ulat)
303       CALL output_field("ulon",f_buf_ulon)
304       CALL output_field("ulat",f_buf_ulat)
305
306       CALL output_field("ps",f_ps)
307       CALL output_field("dps",f_dps)
308       CALL output_field("mass",f_mass)
309       CALL output_field("dmass",f_dmass)
310       CALL output_field("vort",f_qv)
311       CALL output_field("theta",f_theta)
312       CALL output_field("exner",f_pk)
313       CALL output_field("pv",f_qv)
314 
315    END IF
316   
317    !    CALL check_mass_conservation(f_ps,f_dps)
318    CALL trace_end("caldyn")
319!!$OMP BARRIER
320   
321END SUBROUTINE caldyn
322
323SUBROUTINE compute_planetvel(planetvel)
324  USE wind_mod
325  REAL(rstd),INTENT(OUT)  :: planetvel(iim*3*jjm)
326  REAL(rstd) :: ulon(iim*3*jjm)
327  REAL(rstd) :: ulat(iim*3*jjm)
328  REAL(rstd) :: lon,lat
329  INTEGER :: ij
330  DO ij=ij_begin_ext,ij_end_ext
331     CALL xyz2lonlat(xyz_e(ij+u_right,:),lon,lat)
332     ulon(ij+u_right)=a*omega*cos(lat)
333     ulat(ij+u_right)=0
334
335     CALL xyz2lonlat(xyz_e(ij+u_lup,:),lon,lat)
336     ulon(ij+u_lup)=a*omega*cos(lat)
337     ulat(ij+u_lup)=0
338
339     CALL xyz2lonlat(xyz_e(ij+u_ldown,:),lon,lat)
340     ulon(ij+u_ldown)=a*omega*cos(lat)
341     ulat(ij+u_ldown)=0
342  END DO
343  CALL compute_wind2D_perp_from_lonlat_compound(ulon, ulat, planetvel)
344END SUBROUTINE compute_planetvel
345   
346SUBROUTINE compute_pvort(ps,u,theta_rhodz, rhodz,theta,qu,qv)
347  USE icosa
348  USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass
349  USE exner_mod
350  USE trace
351  USE omp_para
352  IMPLICIT NONE
353  REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
354  REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
355  REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm)
356  REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
357  REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm)
358  REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
359  REAL(rstd),INTENT(OUT) :: qv(iim*2*jjm,llm)
360 
361  INTEGER :: i,j,ij,l
362  REAL(rstd) :: etav,hv, m
363!  REAL(rstd) :: qv(2*iim*jjm,llm)     ! potential velocity 
364 
365  CALL trace_start("compute_pvort") 
366
367  IF(caldyn_eta==eta_mass) THEN
368!     CALL wait_message(req_ps) 
369  ELSE
370!     CALL wait_message(req_mass)
371  END IF
372!  CALL wait_message(req_theta_rhodz)
373
374  IF(caldyn_eta==eta_mass) THEN ! Compute mass & theta
375     DO l = ll_begin,ll_end
376!        CALL test_message(req_u)
377!DIR$ SIMD
378        DO ij=ij_begin_ext,ij_end_ext
379           m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
380           rhodz(ij,l) = m
381           theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
382        ENDDO
383     ENDDO
384  ELSE ! Compute only theta
385     DO l = ll_begin,ll_end
386!        CALL test_message(req_u)
387!DIR$ SIMD
388       DO ij=ij_begin_ext,ij_end_ext
389         theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
390       ENDDO
391     ENDDO
392  END IF
393
394!  CALL wait_message(req_u)   
395 
396!!! Compute shallow-water potential vorticity
397  DO l = ll_begin,ll_end
398!DIR$ SIMD
399     DO ij=ij_begin_ext,ij_end_ext
400          etav= 1./Av(ij+z_up)*(  ne_rup        * u(ij+u_rup,l)        * de(ij+u_rup)         &
401                                + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
402                                - ne_lup        * u(ij+u_lup,l)        * de(ij+u_lup) )                               
403
404          hv =  Riv2(ij,vup)          * rhodz(ij,l)            &
405              + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
406              + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
407     
408          qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
409         
410          etav = 1./Av(ij+z_down)*(  ne_ldown         * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
411                                   + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
412                                   - ne_rdown         * u(ij+u_rdown,l)          * de(ij+u_rdown) )
413       
414          hv =  Riv2(ij,vdown)        * rhodz(ij,l)          &
415              + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
416              + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
417                       
418          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
419
420      ENDDO
421
422!DIR$ SIMD
423      DO ij=ij_begin,ij_end
424         qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l)) 
425         qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l)) 
426         qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l)) 
427      END DO
428     
429   ENDDO
430
431   CALL trace_end("compute_pvort")
432                                                       
433  END SUBROUTINE compute_pvort
434 
435  SUBROUTINE compute_geopot(ps,rhodz,theta, pk,geopot) 
436  USE icosa
437  USE disvert_mod
438  USE exner_mod
439  USE trace
440  USE omp_para
441  IMPLICIT NONE
442    REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
443    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
444    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)    ! potential temperature
445    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)       ! Exner function
446    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential
447   
448    INTEGER :: i,j,ij,l
449    REAL(rstd) :: p_ik, exner_ik
450
451    CALL trace_start("compute_geopot")
452
453    IF(caldyn_eta==eta_mass) THEN
454
455!!! Compute exner function and geopotential
456       DO l = 1,llm
457!         !$OMP DO SCHEDULE(STATIC)
458          !DIR$ SIMD
459          DO ij=ij_begin_ext,ij_end_ext         
460             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
461             !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j))
462             exner_ik = cpp * (p_ik/preff) ** kappa
463             pk(ij,l) = exner_ik
464             ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
465             geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik
466          ENDDO
467       ENDDO
468
469    ELSE 
470       ! We are using a Lagrangian vertical coordinate
471       ! Pressure must be computed first top-down (temporarily stored in pk)
472       ! Then Exner pressure and geopotential are computed bottom-up
473       ! Notice that the computation below should work also when caldyn_eta=eta_mass         
474
475       IF(boussinesq) THEN ! compute only geopotential : pressure pk will be computed in compute_caldyn_horiz
476          ! specific volume 1 = dphi/g/rhodz
477          DO l = 1,llm
478!             !$OMP DO SCHEDULE(STATIC)
479             !DIR$ SIMD
480             DO ij=ij_begin_ext,ij_end_ext         
481                geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
482             ENDDO
483          ENDDO
484       ELSE ! non-Boussinesq, compute geopotential and Exner pressure
485          ! uppermost layer
486          !DIR$ SIMD
487          DO ij=ij_begin_ext,ij_end_ext         
488             pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
489          END DO
490          ! other layers
491          DO l = llm-1, 1, -1
492
493!          !$OMP DO SCHEDULE(STATIC)
494             !DIR$ SIMD
495             DO ij=ij_begin_ext,ij_end_ext         
496                pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1))
497             END DO
498          END DO
499          ! surface pressure (for diagnostics)
500          DO ij=ij_begin_ext,ij_end_ext         
501             ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
502          END DO
503
504          ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
505          DO l = 1,llm
506
507!             !$OMP DO SCHEDULE(STATIC)
508             !DIR$ SIMD
509             DO ij=ij_begin_ext,ij_end_ext         
510                p_ik = pk(ij,l)
511                exner_ik = cpp * (p_ik/preff) ** kappa
512                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
513                pk(ij,l) = exner_ik
514             ENDDO
515          ENDDO
516       END IF
517
518    END IF
519
520  CALL trace_end("compute_geopot")
521
522  END SUBROUTINE compute_geopot
523
524  SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,convm, dtheta_rhodz, du)
525  USE icosa
526  USE disvert_mod
527  USE exner_mod
528  USE trace
529  USE omp_para
530  IMPLICIT NONE
531    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)    ! prognostic "velocity"
532    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
533    REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm)
534    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature
535    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function
536    REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential
537
538    REAL(rstd),INTENT(OUT) :: hflux(iim*3*jjm,llm) ! hflux in kg/s
539    REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence
540    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
541    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
542
543    REAL(rstd) :: cor_NT(iim*jjm,llm)  ! NT coriolis force u.(du/dPhi)
544    REAL(rstd) :: urel(3*iim*jjm,llm)  ! relative velocity
545    REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux
546    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
547   
548    INTEGER :: i,j,ij,l
549    REAL(rstd) :: ww,uu
550
551    CALL trace_start("compute_caldyn_horiz")
552
553!    CALL wait_message(req_theta_rhodz)
554
555  DO l = ll_begin, ll_end
556!!!  Compute mass and theta fluxes
557!    IF (caldyn_conserv==energy) CALL test_message(req_qu)
558!DIR$ SIMD
559    DO ij=ij_begin_ext,ij_end_ext         
560        hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
561        hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
562        hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
563
564        Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*hflux(ij+u_right,l)
565        Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l)
566        Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*hflux(ij+u_ldown,l)
567    ENDDO
568
569!!! compute horizontal divergence of fluxes
570!DIR$ SIMD
571    DO ij=ij_begin,ij_end         
572        ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
573        convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
574                                ne_rup*hflux(ij+u_rup,l)       +  & 
575                                ne_lup*hflux(ij+u_lup,l)       +  & 
576                                ne_left*hflux(ij+u_left,l)     +  & 
577                                ne_ldown*hflux(ij+u_ldown,l)   +  & 
578                                ne_rdown*hflux(ij+u_rdown,l))   
579
580        ! signe ? attention d (rho theta dz)
581        ! dtheta_rhodz =  -div(flux.theta)
582        dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right,l)    +  &
583                                 ne_rup*Ftheta(ij+u_rup,l)        +  & 
584                                 ne_lup*Ftheta(ij+u_lup,l)        +  & 
585                                 ne_left*Ftheta(ij+u_left,l)      +  & 
586                                 ne_ldown*Ftheta(ij+u_ldown,l)    +  & 
587                                 ne_rdown*Ftheta(ij+u_rdown,l))
588    ENDDO
589
590 END DO
591
592!!! Compute potential vorticity (Coriolis) contribution to du
593 
594    SELECT CASE(caldyn_conserv)
595    CASE(energy) ! energy-conserving TRiSK
596
597!       CALL wait_message(req_qu)
598
599        DO l=ll_begin,ll_end
600!DIR$ SIMD
601          DO ij=ij_begin,ij_end         
602   
603                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
604                     wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
605                     wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
606                     wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
607                     wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
608                     wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+         &
609                     wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+         &
610                     wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         &
611                     wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+             &
612                     wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))
613                du(ij+u_right,l) = .5*uu/de(ij+u_right)
614               
615                uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
616                     wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
617                     wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
618                     wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
619                     wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
620                     wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) +          &
621                     wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) +              &
622                     wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) +              &
623                     wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) +            &
624                     wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))
625                du(ij+u_lup,l) = .5*uu/de(ij+u_lup)
626
627               
628                uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
629                     wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
630                     wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
631                     wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
632                     wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
633                     wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) +          &
634                     wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) +        &
635                     wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      &
636                     wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) +      &
637                     wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))
638                du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown)
639               
640          ENDDO
641       ENDDO
642       
643    CASE(enstrophy) ! enstrophy-conserving TRiSK
644 
645        DO l=ll_begin,ll_end
646!DIR$ SIMD
647          DO ij=ij_begin,ij_end         
648
649                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
650                     wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
651                     wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
652                     wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
653                     wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
654                     wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
655                     wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
656                     wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
657                     wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
658                     wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
659                du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right)
660               
661               
662                uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
663                     wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
664                     wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
665                     wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
666                     wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
667                     wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
668                     wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
669                     wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
670                     wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
671                     wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
672                du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup)
673
674                uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
675                     wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
676                     wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
677                     wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
678                     wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
679                     wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
680                     wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
681                     wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
682                     wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
683                     wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
684                du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown)
685     
686          ENDDO
687       ENDDO
688       
689    CASE DEFAULT
690       STOP
691    END SELECT
692 
693!!! Compute bernouilli term = Kinetic Energy + geopotential
694    IF(boussinesq) THEN
695       ! first use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure)
696       ! uppermost layer
697       !DIR$ SIMD
698       DO ij=ij_begin_ext,ij_end_ext         
699          pk(ij,llm) = ptop + (.5*g)*theta(ij,llm)*rhodz(ij,llm)
700       END DO
701       ! other layers
702       DO l = llm-1, 1, -1
703!          !$OMP DO SCHEDULE(STATIC)
704          !DIR$ SIMD
705          DO ij=ij_begin_ext,ij_end_ext         
706             pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l)*rhodz(ij,l)+theta(ij,l+1)*rhodz(ij,l+1))
707          END DO
708       END DO
709       ! surface pressure (for diagnostics) FIXME
710       ! DO ij=ij_begin_ext,ij_end_ext         
711       !   ps(ij) = pk(ij,1) + (.5*g)*theta(ij,1)*rhodz(ij,1)
712       ! END DO
713       ! now pk contains the Lagrange multiplier (pressure)
714
715       DO l=ll_begin,ll_end
716!DIR$ SIMD
717          DO ij=ij_begin,ij_end         
718               
719                berni(ij,l) = pk(ij,l) + &
720                       1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
721                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
722                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
723                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
724                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
725                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
726                ! from now on pk contains the vertically-averaged geopotential
727                pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
728          ENDDO
729       ENDDO
730
731    ELSE ! compressible
732
733       DO l=ll_begin,ll_end
734!DIR$ SIMD
735            DO ij=ij_begin,ij_end         
736               
737                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
738                     + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
739                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
740                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
741                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
742                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
743                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
744          ENDDO
745       ENDDO
746
747    END IF ! Boussinesq/compressible
748
749!!! Add gradients of Bernoulli and Exner functions to du
750    DO l=ll_begin,ll_end
751!DIR$ SIMD
752       DO ij=ij_begin,ij_end         
753       
754             du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             &
755                               0.5*(theta(ij,l)+theta(ij+t_right,l))                &
756                                  *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    &
757                                   + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) )
758       
759       
760             du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup) * (                  &
761                               0.5*(theta(ij,l)+theta(ij+t_lup,l))                  &
762                                  *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       &
763                                   + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
764               
765             du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * (             &
766                               0.5*(theta(ij,l)+theta(ij+t_ldown,l))                &
767                                  *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     &
768                                   + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
769
770       ENDDO
771    ENDDO
772 
773    CALL trace_end("compute_caldyn_horiz")
774
775END SUBROUTINE compute_caldyn_horiz
776
777SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
778  USE icosa
779  USE disvert_mod
780  USE exner_mod
781  USE trace
782  USE omp_para
783  IMPLICIT NONE
784    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
785    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)
786    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
787    REAL(rstd),INTENT(INOUT)  :: convm(iim*jjm,llm)  ! mass flux convergence
788
789    REAL(rstd),INTENT(INOUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
790    REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1)
791    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
792    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
793    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
794
795! temporary variable   
796    INTEGER :: i,j,ij,l
797    REAL(rstd) :: p_ik, exner_ik
798
799!    REAL(rstd) :: wwuu(iim*3*jjm,llm+1) ! tmp var, don't know why but gain 30% on the whole code in opemp
800                                        ! need to be understood
801
802!    wwuu=wwuu_out
803  CALL trace_start("compute_caldyn_vert")
804
805!!$OMP BARRIER   
806!!! cumulate mass flux convergence from top to bottom
807  DO  l = llm-1, 1, -1
808!    IF (caldyn_conserv==energy) CALL test_message(req_qu)
809
810!!$OMP DO SCHEDULE(STATIC)
811!DIR$ SIMD
812    DO ij=ij_begin,ij_end         
813        convm(ij,l) = convm(ij,l) + convm(ij,l+1)
814    ENDDO
815  ENDDO
816
817! IMPLICIT FLUSH on convm
818!!!!!!!!!!!!!!!!!!!!!!!!!
819
820! compute dps
821  IF (omp_first) THEN
822!DIR$ SIMD
823     DO ij=ij_begin,ij_end         
824        ! dps/dt = -int(div flux)dz
825        dps(ij) = convm(ij,1) * g 
826    ENDDO
827  ENDIF
828 
829!!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC)
830  DO l=ll_beginp1,ll_end
831!    IF (caldyn_conserv==energy) CALL test_message(req_qu)
832!DIR$ SIMD
833      DO ij=ij_begin,ij_end         
834        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
835        ! => w>0 for upward transport
836        wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l ) 
837    ENDDO
838  ENDDO
839
840  DO l=ll_begin,ll_endm1
841!DIR$ SIMD
842     DO ij=ij_begin,ij_end         
843        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -   0.5 * (  wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1)))
844    ENDDO
845  ENDDO
846 
847  DO l=ll_beginp1,ll_end
848!DIR$ SIMD
849      DO ij=ij_begin,ij_end         
850        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  +   0.5 * ( wflux(ij,l  ) * (theta(ij,l-1) + theta(ij,l) ) )
851    ENDDO
852  ENDDO
853 
854! Compute vertical transport
855  DO l=ll_beginp1,ll_end
856!DIR$ SIMD
857      DO ij=ij_begin,ij_end         
858        wwuu(ij+u_right,l) = 0.5*( wflux(ij,l) + wflux(ij+t_right,l)) * (u(ij+u_right,l) - u(ij+u_right,l-1))
859        wwuu(ij+u_lup,l) = 0.5* ( wflux(ij,l) + wflux(ij+t_lup,l)) * (u(ij+u_lup,l) - u(ij+u_lup,l-1))
860        wwuu(ij+u_ldown,l) = 0.5*( wflux(ij,l) + wflux(ij+t_ldown,l)) * (u(ij+u_ldown,l) - u(ij+u_ldown,l-1))
861     ENDDO
862   ENDDO
863
864 !--> flush wwuu 
865 ! !$OMP BARRIER
866
867! Add vertical transport to du
868  DO l=ll_begin,ll_end
869!DIR$ SIMD
870     DO ij=ij_begin,ij_end         
871        du(ij+u_right, l )   = du(ij+u_right,l)  - (wwuu(ij+u_right,l+1)+ wwuu(ij+u_right,l)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
872        du(ij+u_lup, l )     = du(ij+u_lup,l)    - (wwuu(ij+u_lup,l+1)  + wwuu(ij+u_lup,l))   / (rhodz(ij,l)+rhodz(ij+t_lup,l))
873        du(ij+u_ldown, l )   = du(ij+u_ldown,l)  - (wwuu(ij+u_ldown,l+1)+ wwuu(ij+u_ldown,l)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
874    ENDDO
875  ENDDO     
876
877!  DO l=ll_beginp1,ll_end
878!!DIR$ SIMD
879!      DO ij=ij_begin,ij_end         
880!        wwuu_out(ij+u_right,l) = wwuu(ij+u_right,l)
881!        wwuu_out(ij+u_lup,l)   = wwuu(ij+u_lup,l)
882!        wwuu_out(ij+u_ldown,l) = wwuu(ij+u_ldown,l)
883!     ENDDO
884!   ENDDO
885 
886  CALL trace_end("compute_caldyn_vert")
887
888  END SUBROUTINE compute_caldyn_vert
889
890!-------------------------------- Diagnostics ----------------------------
891
892  SUBROUTINE check_mass_conservation(f_ps,f_dps)
893  USE icosa
894  USE mpipara
895  USE omp_para, ONLY: omp_master
896  IMPLICIT NONE
897    TYPE(t_field),POINTER :: f_ps(:)
898    TYPE(t_field),POINTER :: f_dps(:)
899    REAL(rstd),POINTER :: ps(:)
900    REAL(rstd),POINTER :: dps(:)
901    REAL(rstd) :: mass_tot,dmass_tot
902    INTEGER :: ind,i,j,ij
903   
904    mass_tot=0
905    dmass_tot=0
906   
907    CALL transfert_request(f_dps,req_i1)
908    CALL transfert_request(f_ps,req_i1)
909
910    DO ind=1,ndomain
911      CALL swap_dimensions(ind)
912      CALL swap_geometry(ind)
913
914      ps=f_ps(ind)
915      dps=f_dps(ind)
916
917      DO j=jj_begin,jj_end
918        DO i=ii_begin,ii_end
919          ij=(j-1)*iim+i
920          IF (domain(ind)%own(i,j)) THEN
921            mass_tot=mass_tot+ps(ij)*Ai(ij)/g
922            dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
923          ENDIF
924        ENDDO
925      ENDDO
926   
927    ENDDO
928    IF (is_mpi_root .AND. omp_master) PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot       
929
930  END SUBROUTINE check_mass_conservation 
931 
932  SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
933       f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p)
934    USE icosa
935    USE vorticity_mod
936    USE theta2theta_rhodz_mod
937    USE pression_mod
938    USE omega_mod
939    USE write_field
940    USE vertical_interp_mod
941    USE wind_mod
942    TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_q(:), f_dps(:), &
943         f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:)
944   
945    REAL(rstd) :: out_pression_lev
946    CHARACTER(LEN=255) :: str_pression
947    CHARACTER(LEN=255) :: physics_type
948   
949    out_pression_level=0
950    CALL getin("out_pression_level",out_pression_level) 
951    WRITE(str_pression,*) INT(out_pression_level/100)
952    str_pression=ADJUSTL(str_pression)
953   
954    CALL writefield("ps",f_ps)
955    CALL writefield("dps",f_dps)
956    CALL writefield("phis",f_phis)
957    CALL vorticity(f_u,f_buf_v)
958    CALL writefield("vort",f_buf_v)
959
960    CALL w_omega(f_ps, f_u, f_buf_i)
961    CALL writefield('omega', f_buf_i)
962    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
963      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
964      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
965    ENDIF
966   
967    ! Temperature
968!    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; ! FIXME
969     
970    CALL getin('physics',physics_type)
971    IF (TRIM(physics_type)=='dcmip') THEN
972      CALL Tv2T(f_buf_i,f_q,f_buf1_i) 
973      CALL writefield("T",f_buf1_i)
974      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
975        CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
976        CALL writefield("T"//TRIM(str_pression),f_buf_s)
977      ENDIF
978    ELSE
979      CALL writefield("T",f_buf_i)
980      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
981        CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
982        CALL writefield("T"//TRIM(str_pression),f_buf_s)
983      ENDIF
984    ENDIF
985   
986    ! velocity components
987    CALL un2ulonlat(f_u, f_buf1_i, f_buf2_i)
988    CALL writefield("ulon",f_buf1_i)
989    CALL writefield("ulat",f_buf2_i)
990
991    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
992      CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
993      CALL writefield("ulon"//TRIM(str_pression),f_buf_s)
994      CALL vertical_interp(f_ps,f_buf2_i,f_buf_s,out_pression_level)
995      CALL writefield("ulat"//TRIM(str_pression),f_buf_s)
996    ENDIF
997   
998    ! geopotential ! FIXME
999    CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
1000    CALL writefield("p",f_buf_p)
1001    CALL writefield("phi",f_geopot)   ! geopotential
1002    CALL writefield("theta",f_buf1_i) ! potential temperature
1003    CALL writefield("pk",f_buf2_i)    ! Exner pressure
1004 
1005  END SUBROUTINE write_output_fields
1006 
1007  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi) 
1008    USE field_mod
1009    USE pression_mod
1010    USE exner_mod
1011    USE geopotential_mod
1012    USE theta2theta_rhodz_mod
1013    TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN
1014         f_pks(:), f_p(:), f_theta(:), f_pk(:), f_phi(:)               ! OUT
1015    REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:), theta_rhodz(:,:), &
1016         phi(:,:), phis(:), ps(:), pks(:)
1017    INTEGER :: ind
1018
1019    DO ind=1,ndomain
1020       IF (.NOT. assigned_domain(ind)) CYCLE
1021       CALL swap_dimensions(ind)
1022       CALL swap_geometry(ind)
1023       ps = f_ps(ind)
1024       p  = f_p(ind)
1025       CALL compute_pression(ps,p,0)
1026       pk = f_pk(ind)
1027       pks = f_pks(ind)
1028       CALL compute_exner(ps,p,pks,pk,0)
1029       theta_rhodz = f_theta_rhodz(ind)
1030       theta = f_theta(ind)
1031       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
1032       phis = f_phis(ind)
1033       phi = f_phi(ind)
1034       CALL compute_geopotential(phis,pks,pk,theta,phi,0)
1035    END DO
1036
1037  END SUBROUTINE thetarhodz2geopot
1038 
1039  SUBROUTINE Tv2T(f_Tv, f_q, f_T)
1040  USE icosa
1041  IMPLICIT NONE
1042    TYPE(t_field), POINTER :: f_TV(:)
1043    TYPE(t_field), POINTER :: f_q(:)
1044    TYPE(t_field), POINTER :: f_T(:)
1045   
1046    REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)
1047    INTEGER :: ind
1048   
1049    DO ind=1,ndomain
1050       IF (.NOT. assigned_domain(ind)) CYCLE
1051       CALL swap_dimensions(ind)
1052       CALL swap_geometry(ind)
1053       Tv=f_Tv(ind)
1054       q=f_q(ind)
1055       T=f_T(ind)
1056       T=Tv/(1+0.608*q(:,:,1))
1057    END DO
1058   
1059  END SUBROUTINE Tv2T
1060 
1061END MODULE caldyn_gcm_mod
Note: See TracBrowser for help on using the repository browser.