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

Last change on this file since 802 was 735, checked in by dubos, 6 years ago

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