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

Last change on this file since 171 was 171, checked in by ymipsl, 11 years ago
  • XIOS integration -

Compiling with "-with_xios" option. Adapt path to find XIOS library (arch.path)
Retro-compatible with the old output. If xios is not present, dynamico will use the standard writefield function.
Need to have the iodef.xml configuration file in the exec directory

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