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

Last change on this file since 926 was 920, checked in by dubos, 5 years ago

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