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

Last change on this file since 187 was 186, checked in by ymipsl, 10 years ago

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

YM
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.