source: codes/icosagcm/trunk/src/caldyn_gcm.f90 @ 284

Last change on this file since 284 was 266, checked in by ymipsl, 10 years ago

Synchronize trunk and Saturn branch.
Merge modification from Saturn branch to trunk

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