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

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

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