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

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

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