source: codes/icosagcm/devel/src/dynamics/caldyn_gcm.F90 @ 628

Last change on this file since 628 was 628, checked in by dubos, 7 years ago

devel : new "conservative" variant for vertical transport of momentum
File size: 12.9 KB
Line 
1MODULE caldyn_gcm_mod
2  USE icosa
3  USE transfert_mod
4  USE caldyn_kernels_hevi_mod
5  USE caldyn_kernels_base_mod
6  USE caldyn_kernels_mod
7  USE omp_para
8  USE mpipara
9  IMPLICIT NONE
10  PRIVATE
11
12  PUBLIC init_caldyn, caldyn_BC, caldyn
13 
14CONTAINS
15 
16  SUBROUTINE init_caldyn
17    CHARACTER(len=255) :: def
18    INTEGER            :: ind
19    REAL(rstd),POINTER :: planetvel(:)
20
21    hydrostatic=.TRUE.
22    CALL getin("hydrostatic",hydrostatic)
23 
24    dysl=.NOT.hydrostatic ! dysl defaults to .FALSE. if hydrostatic, .TRUE. if NH                                             
25    CALL getin("dysl",dysl)
26    dysl_geopot=dysl
27    CALL getin("dysl_geopot",dysl_geopot)
28    dysl_caldyn_fast=dysl
29    CALL getin("dysl_caldyn_fast",dysl_caldyn_fast)
30    dysl_slow_hydro=dysl
31    CALL getin("dysl_slow_hydro",dysl_slow_hydro)
32    dysl_pvort_only=dysl
33    CALL getin("dysl_pvort_only",dysl_pvort_only)
34    dysl_caldyn_coriolis=dysl
35    CALL getin("dysl_caldyn_coriolis",dysl_caldyn_coriolis)
36    dysl_caldyn_vert=dysl
37    CALL getin("dysl_caldyn_vert",dysl_caldyn_vert)
38
39    def='advective'
40    CALL getin('caldyn_vert',def)
41    SELECT CASE(TRIM(def))
42    CASE('advective')
43       caldyn_vert_variant=caldyn_vert_noncons
44    CASE('conservative')
45       IF(.NOT. dysl_caldyn_vert) THEN
46          IF (is_mpi_root) PRINT *,'caldyn_vert=cons requires dysl_caldyn_vert=.TRUE.'
47          STOP
48       END IF
49       caldyn_vert_variant=caldyn_vert_cons
50    CASE DEFAULT
51       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_vert : <', &
52            TRIM(def),'> options are <conservative>, <advective>'
53       STOP
54    END SELECT
55    IF (is_master) PRINT *, 'caldyn_vert=',def
56   
57    def='energy'
58    CALL getin('caldyn_conserv',def)
59    SELECT CASE(TRIM(def))
60    CASE('energy')
61       caldyn_conserv=energy
62    CASE('enstrophy')
63       caldyn_conserv=enstrophy
64    CASE DEFAULT
65       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', &
66            TRIM(def),'> options are <energy>, <enstrophy>'
67       STOP
68    END SELECT
69    IF (is_master) PRINT *, 'caldyn_conserv=',def
70
71    nqdyn=1 ! default value
72    physics_thermo = thermo_none
73
74    def='theta'
75    CALL getin('thermo',def)
76    SELECT CASE(TRIM(def))
77    CASE('boussinesq')
78       boussinesq=.TRUE.
79       caldyn_thermo=thermo_boussinesq
80       IF(.NOT. hydrostatic) THEN
81          PRINT *, 'thermo=boussinesq and hydrostatic=.FALSE. : Non-hydrostatic boussinesq equations are not supported'
82          STOP
83       END IF
84    CASE('theta')
85       caldyn_thermo=thermo_theta
86       physics_thermo=thermo_dry
87    CASE('entropy')
88       caldyn_thermo=thermo_entropy
89       physics_thermo=thermo_dry
90    CASE('theta_fake_moist')
91       caldyn_thermo=thermo_theta
92       physics_thermo=thermo_fake_moist
93    CASE('entropy_fake_moist')
94       caldyn_thermo=thermo_entropy
95       physics_thermo=thermo_fake_moist
96    CASE('moist')
97       caldyn_thermo=thermo_moist_debug
98       physics_thermo=thermo_moist
99       nqdyn = 2
100    CASE DEFAULT
101       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_thermo : <', &
102            TRIM(def),'> options are <theta>, <entropy>'
103       STOP
104    END SELECT
105
106    IF(is_master) THEN
107       SELECT CASE(caldyn_thermo)
108       CASE(thermo_theta)
109          PRINT *, 'caldyn_thermo = thermo_theta'
110       CASE(thermo_entropy)
111          PRINT *, 'caldyn_thermo = thermo_entropy'
112       CASE(thermo_moist_debug)
113          PRINT *, 'caldyn_thermo = thermo_moist_debug'
114       CASE DEFAULT
115          STOP
116       END SELECT
117
118       SELECT CASE(physics_thermo)
119       CASE(thermo_dry)
120          PRINT *, 'physics_thermo = thermo_dry'
121       CASE(thermo_fake_moist)
122          PRINT *, 'physics_thermo = thermo_fake_moist'
123       CASE(thermo_moist)
124          PRINT *, 'physics_thermo = thermo_moist'
125       END SELECT
126
127       PRINT *, 'nqdyn =', nqdyn
128    END IF
129
130    CALL allocate_caldyn
131
132    DO ind=1,ndomain
133       IF (.NOT. assigned_domain(ind)) CYCLE
134       CALL swap_dimensions(ind)
135       CALL swap_geometry(ind)
136       planetvel=f_planetvel(ind)
137       CALL compute_planetvel(planetvel)
138    END DO
139
140  END SUBROUTINE init_caldyn
141
142  SUBROUTINE allocate_caldyn
143    CALL allocate_field(f_out_u,field_u,type_real,llm) 
144    CALL allocate_field(f_qu,field_u,type_real,llm) 
145    CALL allocate_field(f_qv,field_z,type_real,llm) 
146    CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk')
147    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu')   
148    CALL allocate_field(f_planetvel,  field_u,type_real, name='planetvel') ! planetary velocity at r=a
149    IF(.NOT.hydrostatic) THEN
150       CALL allocate_field(f_Fel,      field_u,type_real,llm+1,name='F_el')
151       CALL allocate_field(f_gradPhi2, field_t,type_real,llm+1,name='gradPhi2')
152       CALL allocate_field(f_wil,      field_t,type_real,llm+1,name='w_il')
153       CALL allocate_field(f_Wetadot,  field_t,type_real,llm,name='W_etadot')
154    END IF
155  END SUBROUTINE allocate_caldyn
156
157  SUBROUTINE caldyn_BC(f_phis, f_geopot, f_wflux)
158    TYPE(t_field),POINTER :: f_phis(:)
159    TYPE(t_field),POINTER :: f_geopot(:)
160    TYPE(t_field),POINTER :: f_wflux(:)
161    REAL(rstd),POINTER  :: phis(:)
162    REAL(rstd),POINTER  :: wflux(:,:)
163    REAL(rstd),POINTER  :: geopot(:,:)
164    REAL(rstd),POINTER  :: wwuu(:,:)
165
166    INTEGER :: ind,i,j,ij,l
167
168    IF (is_omp_first_level) THEN
169       DO ind=1,ndomain
170          IF (.NOT. assigned_domain(ind)) CYCLE
171          CALL swap_dimensions(ind)
172          CALL swap_geometry(ind)
173          geopot=f_geopot(ind)
174          phis=f_phis(ind)
175          wflux=f_wflux(ind)
176          wwuu=f_wwuu(ind)
177         
178          DO ij=ij_begin_ext,ij_end_ext
179              ! lower BCs : geopot=phis, wflux=0, wwuu=0
180              geopot(ij,1) = phis(ij)
181              wflux(ij,1) = 0.
182              wwuu(ij+u_right,1)=0   
183              wwuu(ij+u_lup,1)=0   
184              wwuu(ij+u_ldown,1)=0
185              ! top BCs : wflux=0, wwuu=0
186              wflux(ij,llm+1)  = 0.
187              wwuu(ij+u_right,llm+1)=0   
188              wwuu(ij+u_lup,llm+1)=0   
189              wwuu(ij+u_ldown,llm+1)=0
190          ENDDO
191       END DO
192    ENDIF
193
194    !$OMP BARRIER
195  END SUBROUTINE caldyn_BC
196   
197  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
198       f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
199    USE observable_mod
200    USE disvert_mod, ONLY : caldyn_eta, eta_mass
201    USE trace
202    LOGICAL,INTENT(IN)    :: write_out
203    TYPE(t_field),POINTER :: f_phis(:)
204    TYPE(t_field),POINTER :: f_ps(:)
205    TYPE(t_field),POINTER :: f_mass(:)
206    TYPE(t_field),POINTER :: f_theta_rhodz(:)
207    TYPE(t_field),POINTER :: f_u(:)
208    TYPE(t_field),POINTER :: f_q(:)
209    TYPE(t_field),POINTER :: f_geopot(:)
210    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
211    TYPE(t_field) :: f_dps(:)
212    TYPE(t_field) :: f_dmass(:)
213    TYPE(t_field) :: f_dtheta_rhodz(:)
214    TYPE(t_field) :: f_du(:)
215   
216    REAL(rstd),POINTER :: ps(:), dps(:)
217    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:)
218    REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
219    REAL(rstd),POINTER :: qu(:,:)
220    REAL(rstd),POINTER :: qv(:,:)
221
222! temporary shared variable
223    REAL(rstd),POINTER  :: theta(:,:,:) 
224    REAL(rstd),POINTER  :: pk(:,:)
225    REAL(rstd),POINTER  :: geopot(:,:)
226    REAL(rstd),POINTER  :: convm(:,:) 
227    REAL(rstd),POINTER  :: wwuu(:,:)
228       
229    INTEGER :: ind
230    LOGICAL,SAVE :: first=.TRUE.
231!$OMP THREADPRIVATE(first)
232   
233    IF (first) THEN
234      first=.FALSE.
235      IF(caldyn_eta==eta_mass) THEN
236         CALL init_message(f_ps,req_i1,req_ps)
237      ELSE
238         CALL init_message(f_mass,req_i1,req_mass)
239      END IF
240      CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
241      CALL init_message(f_u,req_e1_vect,req_u)
242      CALL init_message(f_qu,req_e1_scal,req_qu)
243      ! Overlapping com/compute (deactivated) : MPI messages need to be sent at first call to caldyn
244      ! This is needed only once : the next ones will be sent by timeloop
245!      IF(caldyn_eta==eta_mass) THEN
246!         CALL send_message(f_ps,req_ps)
247!         CALL wait_message(req_ps) 
248!      ELSE
249!         CALL send_message(f_mass,req_mass)
250!         CALL wait_message(req_mass) 
251!      END IF
252    ENDIF
253   
254    CALL trace_start("caldyn")
255
256    IF(caldyn_eta==eta_mass) THEN
257       CALL send_message(f_ps,req_ps) ! COM00
258       CALL wait_message(req_ps) ! COM00
259    ELSE
260       CALL send_message(f_mass,req_mass) ! COM00
261       CALL wait_message(req_mass) ! COM00
262    END IF
263   
264    CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01
265    CALL wait_message(req_theta_rhodz) ! COM01 Moved from caldyn_pvort
266    CALL send_message(f_u,req_u) ! COM02
267    CALL wait_message(req_u) ! COM02
268
269    IF(.NOT.hydrostatic) THEN
270       STOP 'caldyn_gcm may not be used yet when non-hydrostatic'
271    END IF
272
273    SELECT CASE(caldyn_conserv)
274    CASE(energy) ! energy-conserving
275       DO ind=1,ndomain
276          IF (.NOT. assigned_domain(ind)) CYCLE
277          CALL swap_dimensions(ind)
278          CALL swap_geometry(ind)
279          ps=f_ps(ind)
280          u=f_u(ind)
281          theta_rhodz = f_theta_rhodz(ind)
282          mass=f_mass(ind)
283          theta = f_theta(ind)
284          qu=f_qu(ind)
285          qv=f_qv(ind)
286          pk = f_pk(ind)
287          geopot = f_geopot(ind) 
288          hflux=f_hflux(ind)
289          convm = f_dmass(ind)
290          dtheta_rhodz=f_dtheta_rhodz(ind)
291          du=f_du(ind)
292          CALL compute_pvort(ps,u,theta_rhodz(:,:,1), mass,theta,qu,qv) ! COM00 COM01 COM02
293!          CALL compute_theta(ps,theta_rhodz, mass,theta)
294!          CALL compute_pvort_only(u,mass,qu,qv)
295
296          CALL compute_geopot(mass,theta, ps,pk,geopot)
297!          du(:,:)=0.
298!          CALL compute_caldyn_fast(0.,u,mass,theta,pk,geopot,du)
299       ENDDO       
300
301       CALL send_message(f_u,req_u) ! COM02
302       CALL wait_message(req_u) ! COM02
303       CALL send_message(f_qu,req_qu) ! COM03
304       CALL wait_message(req_qu) ! COM03
305
306       DO ind=1,ndomain
307          IF (.NOT. assigned_domain(ind)) CYCLE
308          CALL swap_dimensions(ind)
309          CALL swap_geometry(ind)
310          ps=f_ps(ind)
311          u=f_u(ind)
312          theta_rhodz = f_theta_rhodz(ind)
313          mass=f_mass(ind)
314          theta = f_theta(ind)
315          qu=f_qu(ind)
316          qv=f_qv(ind)
317          pk = f_pk(ind)
318          geopot = f_geopot(ind) 
319          hflux=f_hflux(ind)
320          convm = f_dmass(ind)
321          dtheta_rhodz=f_dtheta_rhodz(ind)
322          du=f_du(ind)
323
324          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz(:,:,1),du)
325!          CALL compute_caldyn_slow_hydro(u,mass,hflux,du, .FALSE.) ! FIXME
326!          CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
327          IF(caldyn_eta==eta_mass) THEN
328             wflux=f_wflux(ind)
329             wwuu=f_wwuu(ind)
330             dps=f_dps(ind)
331             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz(:,:,1), du)
332          END IF
333       ENDDO       
334   
335    CASE(enstrophy) ! enstrophy-conserving
336       DO ind=1,ndomain
337          IF (.NOT. assigned_domain(ind)) CYCLE
338          CALL swap_dimensions(ind)
339          CALL swap_geometry(ind)
340          ps=f_ps(ind)
341          u=f_u(ind)
342          theta_rhodz=f_theta_rhodz(ind)
343          mass=f_mass(ind)
344          theta = f_theta(ind)
345          qu=f_qu(ind)
346          qv=f_qv(ind)
347          CALL compute_pvort(ps,u,theta_rhodz(:,:,1), mass,theta,qu,qv)
348          pk = f_pk(ind)
349          geopot = f_geopot(ind) 
350          CALL compute_geopot(ps,mass,theta, pk,geopot)
351          hflux=f_hflux(ind)
352          convm = f_dmass(ind)
353          dtheta_rhodz=f_dtheta_rhodz(ind)
354          du=f_du(ind)
355          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz(:,:,1),du)
356          IF(caldyn_eta==eta_mass) THEN
357             wflux=f_wflux(ind)
358             wwuu=f_wwuu(ind)
359             dps=f_dps(ind)
360             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
361          END IF
362       ENDDO
363       
364    CASE DEFAULT
365       STOP
366    END SELECT
367
368!$OMP BARRIER
369    !    CALL check_mass_conservation(f_ps,f_dps)
370    CALL trace_end("caldyn")
371!!$OMP BARRIER
372   
373END SUBROUTINE caldyn
374
375!-------------------------------- Diagnostics ----------------------------
376
377  SUBROUTINE check_mass_conservation(f_ps,f_dps)
378    TYPE(t_field),POINTER :: f_ps(:)
379    TYPE(t_field),POINTER :: f_dps(:)
380    REAL(rstd),POINTER :: ps(:)
381    REAL(rstd),POINTER :: dps(:)
382    REAL(rstd) :: mass_tot,dmass_tot
383    INTEGER :: ind,i,j,ij
384   
385    mass_tot=0
386    dmass_tot=0
387   
388    CALL transfert_request(f_dps,req_i1)
389    CALL transfert_request(f_ps,req_i1)
390
391    DO ind=1,ndomain
392      CALL swap_dimensions(ind)
393      CALL swap_geometry(ind)
394
395      ps=f_ps(ind)
396      dps=f_dps(ind)
397
398      DO j=jj_begin,jj_end
399        DO i=ii_begin,ii_end
400          ij=(j-1)*iim+i
401          IF (domain(ind)%own(i,j)) THEN
402            mass_tot=mass_tot+ps(ij)*Ai(ij)/g
403            dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
404          ENDIF
405        ENDDO
406      ENDDO
407   
408    ENDDO
409    IF (is_mpi_root) PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot       
410
411  END SUBROUTINE check_mass_conservation 
412 
413END MODULE caldyn_gcm_mod
Note: See TracBrowser for help on using the repository browser.