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

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

Multi-layer Saint-Venant / Ripa equations - tested with Williamson91.5 & JW06
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    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, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass
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 = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
321              rhodz(ij,l) = m
322              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
323           ENDDO
324        ENDDO
325     ENDDO
326  ELSE ! Compute only theta
327     DO l = ll_begin,ll_end
328        CALL test_message(req_u) 
329        DO j=jj_begin-1,jj_end+1
330           DO i=ii_begin-1,ii_end+1
331              ij=(j-1)*iim+i
332              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
333           ENDDO
334        ENDDO
335     ENDDO
336  END IF
337
338  CALL wait_message(req_u)   
339 
340!!! Compute shallow-water potential vorticity
341  DO l = ll_begin,ll_end
342
343    DO j=jj_begin-1,jj_end+1
344       DO i=ii_begin-1,ii_end+1
345          ij=(j-1)*iim+i
346         
347          etav= 1./Av(ij+z_up)*(  ne_rup        * u(ij+u_rup,l)        * de(ij+u_rup)         &
348                                + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
349                                - ne_lup        * u(ij+u_lup,l)        * de(ij+u_lup) )                               
350
351          hv =  Riv2(ij,vup)          * rhodz(ij,l)            &
352              + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
353              + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
354     
355          qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
356         
357          etav = 1./Av(ij+z_down)*(  ne_ldown         * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
358                                   + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
359                                   - ne_rdown         * u(ij+u_rdown,l)          * de(ij+u_rdown) )
360       
361          hv =  Riv2(ij,vdown)        * rhodz(ij,l)          &
362              + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
363              + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
364                       
365          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
366
367        ENDDO
368      ENDDO
369
370      DO j=jj_begin,jj_end
371         DO i=ii_begin,ii_end
372            ij=(j-1)*iim+i
373            qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l)) 
374            qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l)) 
375            qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l)) 
376         END DO
377      END DO
378     
379   ENDDO
380
381   CALL trace_end("compute_pvort")
382                                                       
383  END SUBROUTINE compute_pvort
384 
385  SUBROUTINE compute_geopot(ps,rhodz,theta,pk,geopot) 
386  USE icosa
387  USE disvert_mod
388  USE exner_mod
389  USE trace
390  USE omp_para
391  IMPLICIT NONE
392    REAL(rstd),INTENT(INOUT) :: ps(iim*jjm)
393    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
394    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)    ! potential temperature
395    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)       ! Exner function
396    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) ! geopotential
397   
398    INTEGER :: i,j,ij,l
399    REAL(rstd) :: p_ik, exner_ik
400
401    CALL trace_start("compute_geopot")
402
403    IF(caldyn_eta==eta_mass) THEN
404   
405!!! Compute exner function and geopotential
406       DO l = 1,llm
407          !$OMP DO SCHEDULE(STATIC)
408          DO j=jj_begin-1,jj_end+1
409             DO i=ii_begin-1,ii_end+1
410                ij=(j-1)*iim+i
411                p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
412                !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j))
413                exner_ik = cpp * (p_ik/preff) ** kappa
414                pk(ij,l) = exner_ik
415                ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
416                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik
417             ENDDO
418          ENDDO
419       ENDDO
420
421    ELSE 
422       ! We are using a Lagrangian vertical coordinate
423       ! Pressure must be computed first top-down (temporarily stored in pk)
424       ! Then Exner pressure and geopotential are computed bottom-up
425       ! Notice that the computation below should work also when caldyn_eta=eta_mass         
426
427       ! uppermost layer
428       DO j=jj_begin-1,jj_end+1
429          DO i=ii_begin-1,ii_end+1
430             ij=(j-1)*iim+i
431             pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
432          END DO
433       END DO
434       ! other layers
435       DO l = llm-1, 1, -1
436          !$OMP DO SCHEDULE(STATIC)
437          DO j=jj_begin-1,jj_end+1
438             DO i=ii_begin-1,ii_end+1
439                ij=(j-1)*iim+i
440                pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1))
441             END DO
442          END DO
443       END DO
444       ! surface pressure (for diagnostics)
445       DO j=jj_begin-1,jj_end+1
446          DO i=ii_begin-1,ii_end+1
447             ij=(j-1)*iim+i
448             ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
449          END DO
450       END DO
451
452       IF(boussinesq) THEN ! specific volume 1 = dphi/g/rhodz
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                   geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
459                ENDDO
460             ENDDO
461          ENDDO
462       ELSE ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
463          DO l = 1,llm
464             !$OMP DO SCHEDULE(STATIC)
465             DO j=jj_begin-1,jj_end+1
466                DO i=ii_begin-1,ii_end+1
467                   ij=(j-1)*iim+i
468                   p_ik = pk(ij,l)
469                   exner_ik = cpp * (p_ik/preff) ** kappa
470                   geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
471                   pk(ij,l) = exner_ik
472                ENDDO
473             ENDDO
474          ENDDO
475       END IF
476
477    END IF
478
479  CALL trace_end("compute_geopot")
480
481  END SUBROUTINE compute_geopot
482
483  SUBROUTINE compute_caldyn_horiz(u,rhodz,qu,theta,pk,geopot, hflux,convm, dtheta_rhodz, du)
484  USE icosa
485  USE disvert_mod
486  USE exner_mod
487  USE trace
488  USE omp_para
489  IMPLICIT NONE
490    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
491    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
492    REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm)
493    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)  ! potential temperature
494    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm) ! Exner function
495    REAL(rstd),INTENT(IN)  :: geopot(iim*jjm,llm+1)    ! geopotential
496
497    REAL(rstd),INTENT(OUT) :: hflux(iim*3*jjm,llm) ! hflux in kg/s
498    REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence
499    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
500    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
501
502    REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux
503    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
504   
505    INTEGER :: i,j,ij,l
506    REAL(rstd) :: ww,uu
507
508    CALL trace_start("compute_caldyn_horiz")
509
510    CALL wait_message(req_theta_rhodz) 
511
512  DO l = ll_begin, ll_end
513!!!  Compute mass and theta fluxes
514    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
515    DO j=jj_begin-1,jj_end+1
516      DO i=ii_begin-1,ii_end+1
517        ij=(j-1)*iim+i
518        hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
519        hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
520        hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
521
522        Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*hflux(ij+u_right,l)
523        Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l)
524        Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*hflux(ij+u_ldown,l)
525      ENDDO
526    ENDDO
527
528!!! compute horizontal divergence of fluxes
529    DO j=jj_begin,jj_end
530      DO i=ii_begin,ii_end
531        ij=(j-1)*iim+i 
532        ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
533        convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
534                                ne_rup*hflux(ij+u_rup,l)       +  & 
535                                ne_lup*hflux(ij+u_lup,l)       +  & 
536                                ne_left*hflux(ij+u_left,l)     +  & 
537                                ne_ldown*hflux(ij+u_ldown,l)   +  & 
538                                ne_rdown*hflux(ij+u_rdown,l))   
539
540        ! signe ? attention d (rho theta dz)
541        ! dtheta_rhodz =  -div(flux.theta)
542        dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right,l)    +  &
543                                 ne_rup*Ftheta(ij+u_rup,l)        +  & 
544                                 ne_lup*Ftheta(ij+u_lup,l)        +  & 
545                                 ne_left*Ftheta(ij+u_left,l)      +  & 
546                                 ne_ldown*Ftheta(ij+u_ldown,l)    +  & 
547                                 ne_rdown*Ftheta(ij+u_rdown,l))
548      ENDDO
549    ENDDO
550
551 END DO
552
553!!! Compute potential vorticity (Coriolis) contribution to du
554 
555    SELECT CASE(caldyn_conserv)
556    CASE(energy) ! energy-conserving TRiSK
557
558       CALL wait_message(req_qu)
559
560        DO l=ll_begin,ll_end
561          DO j=jj_begin,jj_end
562             DO i=ii_begin,ii_end
563                ij=(j-1)*iim+i
564
565                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
566                     wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
567                     wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
568                     wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
569                     wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
570                     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))+         &
571                     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))+         &
572                     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))+         &
573                     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))+             &
574                     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))
575                du(ij+u_right,l) = .5*uu/de(ij+u_right)
576               
577                uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
578                     wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
579                     wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
580                     wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
581                     wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
582                     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)) +          &
583                     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)) +              &
584                     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)) +              &
585                     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)) +            &
586                     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))
587                du(ij+u_lup,l) = .5*uu/de(ij+u_lup)
588
589               
590                uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
591                     wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
592                     wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
593                     wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
594                     wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
595                     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)) +          &
596                     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)) +        &
597                     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)) +      &
598                     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)) +      &
599                     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))
600                du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown)
601               
602             ENDDO
603          ENDDO
604       ENDDO
605       
606    CASE(enstrophy) ! enstrophy-conserving TRiSK
607 
608        DO l=ll_begin,ll_end
609          DO j=jj_begin,jj_end
610             DO i=ii_begin,ii_end
611                ij=(j-1)*iim+i
612
613                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
614                     wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
615                     wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
616                     wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
617                     wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
618                     wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
619                     wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
620                     wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
621                     wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
622                     wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
623                du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right)
624               
625               
626                uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
627                     wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
628                     wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
629                     wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
630                     wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
631                     wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
632                     wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
633                     wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
634                     wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
635                     wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
636                du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup)
637
638                uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
639                     wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
640                     wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
641                     wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
642                     wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
643                     wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
644                     wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
645                     wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
646                     wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
647                     wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
648                du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown)
649     
650             ENDDO
651          ENDDO
652       ENDDO
653       
654    CASE DEFAULT
655       STOP
656    END SELECT
657 
658!!! Compute bernouilli term = Kinetic Energy + geopotential
659    IF(boussinesq) THEN
660
661       DO l=ll_begin,ll_end
662          DO j=jj_begin,jj_end
663             DO i=ii_begin,ii_end
664                ij=(j-1)*iim+i
665               
666                berni(ij,l) = pk(ij,l) + &
667                     + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
668                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
669                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
670                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
671                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
672                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
673                pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
674             ENDDO
675          ENDDO
676       ENDDO
677
678    ELSE
679
680       DO l=ll_begin,ll_end
681          DO j=jj_begin,jj_end
682             DO i=ii_begin,ii_end
683                ij=(j-1)*iim+i
684               
685                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
686                     + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
687                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
688                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
689                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
690                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
691                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
692       
693             ENDDO
694          ENDDO
695       ENDDO
696
697    END IF ! Boussinesq/compressible
698
699!!! Add gradients of Bernoulli and Exner functions to du
700    DO l=ll_begin,ll_end
701       DO j=jj_begin,jj_end
702          DO i=ii_begin,ii_end
703             ij=(j-1)*iim+i
704       
705             du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             &
706                               0.5*(theta(ij,l)+theta(ij+t_right,l))                &
707                                  *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    &
708                                   + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) )
709       
710       
711             du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup) * (                  &
712                               0.5*(theta(ij,l)+theta(ij+t_lup,l))                  &
713                                  *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       &
714                                   + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
715               
716             du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * (             &
717                               0.5*(theta(ij,l)+theta(ij+t_ldown,l))                &
718                                  *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     &
719                                   + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
720
721          ENDDO
722       ENDDO
723    ENDDO
724 
725    CALL trace_end("compute_caldyn_horiz")
726
727END SUBROUTINE compute_caldyn_horiz
728
729SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
730  USE icosa
731  USE disvert_mod
732  USE exner_mod
733  USE trace
734  USE omp_para
735  IMPLICIT NONE
736    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
737    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)
738    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
739    REAL(rstd),INTENT(INOUT)  :: convm(iim*jjm,llm)  ! mass flux convergence
740
741    REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
742    REAL(rstd),INTENT(OUT) :: wwuu(iim*3*jjm,llm+1)
743    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm)
744    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
745    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
746
747! temporary variable   
748    INTEGER :: i,j,ij,l
749    REAL(rstd) :: p_ik, exner_ik
750
751  CALL trace_start("compute_caldyn_vert")
752
753!$OMP BARRIER   
754!!! cumulate mass flux convergence from top to bottom
755  DO  l = llm-1, 1, -1
756    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
757!$OMP DO SCHEDULE(STATIC)
758    DO j=jj_begin,jj_end
759      DO i=ii_begin,ii_end
760        ij=(j-1)*iim+i
761        convm(ij,l) = convm(ij,l) + convm(ij,l+1)
762      ENDDO
763    ENDDO
764  ENDDO
765
766! IMPLICIT FLUSH on convm
767!!!!!!!!!!!!!!!!!!!!!!!!!
768
769! compute dps
770  IF (omp_first) THEN
771    DO j=jj_begin,jj_end
772      DO i=ii_begin,ii_end
773        ij=(j-1)*iim+i
774        ! dps/dt = -int(div flux)dz
775        dps(ij) = convm(ij,1) * g 
776      ENDDO
777    ENDDO
778  ENDIF
779 
780!!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC)
781  DO l=ll_beginp1,ll_end
782    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
783    DO j=jj_begin,jj_end
784      DO i=ii_begin,ii_end
785        ij=(j-1)*iim+i
786        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
787        ! => w>0 for upward transport
788        wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l ) 
789      ENDDO
790    ENDDO
791  ENDDO
792
793  DO l=ll_begin,ll_endm1
794    DO j=jj_begin,jj_end
795      DO i=ii_begin,ii_end
796        ij=(j-1)*iim+i
797        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -   0.5 * (  wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1)))
798      ENDDO
799    ENDDO
800  ENDDO
801 
802  DO l=ll_beginp1,ll_end
803    DO j=jj_begin,jj_end
804      DO i=ii_begin,ii_end
805        ij=(j-1)*iim+i
806        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  +   0.5 * ( wflux(ij,l  ) * (theta(ij,l-1) + theta(ij,l) ) )
807      ENDDO
808    ENDDO
809  ENDDO
810 
811! Compute vertical transport
812  DO l=ll_beginp1,ll_end
813    DO j=jj_begin,jj_end
814      DO i=ii_begin,ii_end
815        ij=(j-1)*iim+i
816        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))
817        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))
818        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))
819       ENDDO
820     ENDDO
821   ENDDO
822
823 !--> flush wwuu 
824 !$OMP BARRIER
825
826! Add vertical transport to du
827  DO l=ll_begin,ll_end
828    DO j=jj_begin,jj_end
829      DO i=ii_begin,ii_end
830        ij=(j-1)*iim+i
831        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))
832        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))
833        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))
834      ENDDO
835    ENDDO
836  ENDDO     
837 
838  CALL trace_end("compute_caldyn_vert")
839
840  END SUBROUTINE compute_caldyn_vert
841
842!-------------------------------- Diagnostics ----------------------------
843
844  SUBROUTINE check_mass_conservation(f_ps,f_dps)
845  USE icosa
846  USE mpipara
847  IMPLICIT NONE
848    TYPE(t_field),POINTER :: f_ps(:)
849    TYPE(t_field),POINTER :: f_dps(:)
850    REAL(rstd),POINTER :: ps(:)
851    REAL(rstd),POINTER :: dps(:)
852    REAL(rstd) :: mass_tot,dmass_tot
853    INTEGER :: ind,i,j,ij
854   
855    mass_tot=0
856    dmass_tot=0
857   
858    CALL transfert_request(f_dps,req_i1)
859    CALL transfert_request(f_ps,req_i1)
860
861    DO ind=1,ndomain
862      CALL swap_dimensions(ind)
863      CALL swap_geometry(ind)
864
865      ps=f_ps(ind)
866      dps=f_dps(ind)
867
868      DO j=jj_begin,jj_end
869        DO i=ii_begin,ii_end
870          ij=(j-1)*iim+i
871          IF (domain(ind)%own(i,j)) THEN
872            mass_tot=mass_tot+ps(ij)*Ai(ij)/g
873            dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
874          ENDIF
875        ENDDO
876      ENDDO
877   
878    ENDDO
879    IF (is_mpi_root) PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot       
880
881  END SUBROUTINE check_mass_conservation 
882 
883  SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
884       f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p)
885    USE icosa
886    USE vorticity_mod
887    USE theta2theta_rhodz_mod
888    USE pression_mod
889    USE omega_mod
890    USE write_field
891    USE vertical_interp_mod
892    USE wind_mod
893    TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_q(:), f_dps(:), &
894         f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:)
895   
896    REAL(rstd) :: out_pression_lev
897    CHARACTER(LEN=255) :: str_pression
898    CHARACTER(LEN=255) :: physics_type
899   
900    out_pression_level=0
901    CALL getin("out_pression_level",out_pression_level) 
902    WRITE(str_pression,*) INT(out_pression_level/100)
903    str_pression=ADJUSTL(str_pression)
904   
905    CALL writefield("ps",f_ps)
906    CALL writefield("dps",f_dps)
907    CALL writefield("phis",f_phis)
908    CALL vorticity(f_u,f_buf_v)
909    CALL writefield("vort",f_buf_v)
910
911    CALL w_omega(f_ps, f_u, f_buf_i)
912    CALL writefield('omega', f_buf_i)
913    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
914      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
915      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
916    ENDIF
917   
918    ! Temperature
919!    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; ! FIXME
920     
921    CALL getin('physics',physics_type)
922    IF (TRIM(physics_type)=='dcmip') THEN
923      CALL Tv2T(f_buf_i,f_q,f_buf1_i) 
924      CALL writefield("T",f_buf1_i)
925      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
926        CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
927        CALL writefield("T"//TRIM(str_pression),f_buf_s)
928      ENDIF
929    ELSE
930      CALL writefield("T",f_buf_i)
931      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
932        CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
933        CALL writefield("T"//TRIM(str_pression),f_buf_s)
934      ENDIF
935    ENDIF
936   
937    ! velocity components
938    CALL un2ulonlat(f_u, f_buf1_i, f_buf2_i)
939    CALL writefield("ulon",f_buf1_i)
940    CALL writefield("ulat",f_buf2_i)
941
942    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
943      CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
944      CALL writefield("ulon"//TRIM(str_pression),f_buf_s)
945      CALL vertical_interp(f_ps,f_buf2_i,f_buf_s,out_pression_level)
946      CALL writefield("ulat"//TRIM(str_pression),f_buf_s)
947    ENDIF
948   
949    ! geopotential ! FIXME
950    CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
951    CALL writefield("p",f_buf_p)
952    CALL writefield("phi",f_geopot)   ! geopotential
953    CALL writefield("theta",f_buf1_i) ! potential temperature
954    CALL writefield("pk",f_buf2_i)    ! Exner pressure
955 
956  END SUBROUTINE write_output_fields
957 
958  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi) 
959    USE field_mod
960    USE pression_mod
961    USE exner_mod
962    USE geopotential_mod
963    USE theta2theta_rhodz_mod
964    TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN
965         f_pks(:), f_p(:), f_theta(:), f_pk(:), f_phi(:)               ! OUT
966    REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:), theta_rhodz(:,:), &
967         phi(:,:), phis(:), ps(:), pks(:)
968    INTEGER :: ind
969
970    DO ind=1,ndomain
971       CALL swap_dimensions(ind)
972       CALL swap_geometry(ind)
973       ps = f_ps(ind)
974       p  = f_p(ind)
975       CALL compute_pression(ps,p,0)
976       pk = f_pk(ind)
977       pks = f_pks(ind)
978       CALL compute_exner(ps,p,pks,pk,0)
979       theta_rhodz = f_theta_rhodz(ind)
980       theta = f_theta(ind)
981       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
982       phis = f_phis(ind)
983       phi = f_phi(ind)
984       CALL compute_geopotential(phis,pks,pk,theta,phi,0)
985    END DO
986
987  END SUBROUTINE thetarhodz2geopot
988 
989  SUBROUTINE Tv2T(f_Tv, f_q, f_T)
990  USE icosa
991  IMPLICIT NONE
992    TYPE(t_field), POINTER :: f_TV(:)
993    TYPE(t_field), POINTER :: f_q(:)
994    TYPE(t_field), POINTER :: f_T(:)
995   
996    REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)
997    INTEGER :: ind
998   
999    DO ind=1,ndomain
1000       CALL swap_dimensions(ind)
1001       CALL swap_geometry(ind)
1002       Tv=f_Tv(ind)
1003       q=f_q(ind)
1004       T=f_T(ind)
1005       T=Tv/(1+0.608*q(:,:,1))
1006    END DO
1007   
1008  END SUBROUTINE Tv2T
1009 
1010END MODULE caldyn_gcm_mod
Note: See TracBrowser for help on using the repository browser.