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

Last change on this file since 162 was 162, checked in by dubos, 11 years ago

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