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

Last change on this file since 354 was 354, checked in by dubos, 9 years ago

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