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

Last change on this file since 190 was 190, checked in by dubos, 10 years ago

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