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

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

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