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

Last change on this file since 327 was 327, checked in by ymipsl, 9 years ago

Merge recent developments from saturn branch onto trunk.

  • lmdz generic physics interface
  • performance improvment on mix mpi/openmp
  • asynchrone and overlaping communication
  • best domain distribution between process and threads
  • ....

This version is compatible with the actual saturn version and the both branches are considered merged on dynamico component.

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