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

Last change on this file since 905 was 885, checked in by dubos, 5 years ago

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