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

Last change on this file since 176 was 174, checked in by ymipsl, 11 years ago

Transform 2 loops on i and j in one loop ij for efficiency, vectorization and future GPU programing

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