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

Last change on this file since 1037 was 928, checked in by dubos, 5 years ago

devel : separate modules for caldyn_vert and caldyn_vert_NH
File size: 14.5 KB
RevLine 
[12]1MODULE caldyn_gcm_mod
[19]2  USE icosa
[151]3  USE transfert_mod
[731]4  USE caldyn_vars_mod
[928]5  USE compute_caldyn_vert_mod, ONLY : compute_caldyn_vert => compute_caldyn_vert_manual
[361]6  USE caldyn_kernels_mod
[586]7  USE omp_para
8  USE mpipara
[350]9  IMPLICIT NONE
[132]10  PRIVATE
11
[362]12  PUBLIC init_caldyn, caldyn_BC, caldyn
13 
[12]14CONTAINS
[15]15 
[98]16  SUBROUTINE init_caldyn
[827]17
[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
[735]63    CASE('energy_gassmann')
64       caldyn_conserv=conserv_gassmann
[126]65    CASE('enstrophy')
[733]66       caldyn_conserv=conserv_enstrophy
[125]67    CASE DEFAULT
[159]68       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', &
[735]69            TRIM(def),'> options are <energy>, <energy_gassmann>, <enstrophy>'
[125]70       STOP
71    END SELECT
[295]72    IF (is_master) PRINT *, 'caldyn_conserv=',def
[125]73
[733]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
[413]88    nqdyn=1 ! default value
[534]89    physics_thermo = thermo_none
[413]90
[401]91    def='theta'
[413]92    CALL getin('thermo',def)
[401]93    SELECT CASE(TRIM(def))
[534]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
[401]101    CASE('theta')
102       caldyn_thermo=thermo_theta
[406]103       physics_thermo=thermo_dry
[401]104    CASE('entropy')
105       caldyn_thermo=thermo_entropy
[406]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')
[413]114       caldyn_thermo=thermo_moist_debug
[406]115       physics_thermo=thermo_moist
[413]116       nqdyn = 2
[401]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
[413]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
[17]147    CALL allocate_caldyn
[861]148   
149    SELECT CASE(grid_type)
150    CASE(grid_ico)
[846]151       ! NB : unused now, would be useful for non-traditional / deep-atmosphere
152       DO ind=1,ndomain
153          IF (.NOT. assigned_domain(ind)) CYCLE
154          CALL swap_dimensions(ind)
155          CALL swap_geometry(ind)
156          planetvel=f_planetvel(ind)
157          CALL compute_planetvel(planetvel)
158       END DO
[861]159    CASE(grid_unst)
160    END SELECT
[179]161
[15]162  END SUBROUTINE init_caldyn
[179]163
[12]164  SUBROUTINE allocate_caldyn
[735]165    CALL allocate_field(f_qu,field_u,type_real,llm, name='qu')
166    CALL allocate_field(f_qv,field_z,type_real,llm, name='qv')
167    CALL allocate_field(f_Kv,field_z,type_real,llm, name='Kv')
168    CALL allocate_field(f_hv,field_z,type_real,llm, name='hv')
[159]169    CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk')
[735]170    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu')
[179]171    CALL allocate_field(f_planetvel,  field_u,type_real, name='planetvel') ! planetary velocity at r=a
[928]172    CALL allocate_field(f_Ftheta,field_u,type_real,llm, name='Ftheta')
173    CALL allocate_field(f_berni, field_t,type_real,llm, name='berni')
[558]174    IF(.NOT.hydrostatic) THEN
175       CALL allocate_field(f_Fel,      field_u,type_real,llm+1,name='F_el')
176       CALL allocate_field(f_gradPhi2, field_t,type_real,llm+1,name='gradPhi2')
177       CALL allocate_field(f_wil,      field_t,type_real,llm+1,name='w_il')
178       CALL allocate_field(f_Wetadot,  field_t,type_real,llm,name='W_etadot')
179    END IF
[12]180  END SUBROUTINE allocate_caldyn
[157]181
[350]182  SUBROUTINE caldyn_BC(f_phis, f_geopot, f_wflux)
[157]183    TYPE(t_field),POINTER :: f_phis(:)
[350]184    TYPE(t_field),POINTER :: f_geopot(:)
[157]185    TYPE(t_field),POINTER :: f_wflux(:)
186    REAL(rstd),POINTER  :: phis(:)
187    REAL(rstd),POINTER  :: wflux(:,:)
188    REAL(rstd),POINTER  :: geopot(:,:)
189    REAL(rstd),POINTER  :: wwuu(:,:)
190
191    INTEGER :: ind,i,j,ij,l
192
[295]193    IF (is_omp_first_level) THEN
[157]194       DO ind=1,ndomain
[186]195          IF (.NOT. assigned_domain(ind)) CYCLE
[157]196          CALL swap_dimensions(ind)
197          CALL swap_geometry(ind)
198          geopot=f_geopot(ind)
199          phis=f_phis(ind)
200          wflux=f_wflux(ind)
201          wwuu=f_wwuu(ind)
202         
[885]203          SELECT CASE(grid_type)
204          CASE(grid_unst)
205             ! FIXME OMP
206             DO ij=1, primal_num
207                geopot(1,ij) = phis(ij)
208                wflux(1,ij)=0.
209                wflux(llm+1,ij)=0.
210             END DO
211             DO ij=1, edge_num
212                wwuu(1,ij)=0.
213                wwuu(llm+1,ij)=0.
214             END DO
215          CASE(grid_ico)
216             DO ij=ij_begin_ext,ij_end_ext
217                ! lower BCs : geopot=phis, wflux=0, wwuu=0
218                geopot(ij,1) = phis(ij)
219                wflux(ij,1) = 0.
220                wwuu(ij+u_right,1)=0   
221                wwuu(ij+u_lup,1)=0   
222                wwuu(ij+u_ldown,1)=0
223                ! top BCs : wflux=0, wwuu=0
224                wflux(ij,llm+1)  = 0.
225                wwuu(ij+u_right,llm+1)=0   
226                wwuu(ij+u_lup,llm+1)=0   
227                wwuu(ij+u_ldown,llm+1)=0
228             END DO
229          END SELECT
[157]230       END DO
231    ENDIF
232
[295]233    !$OMP BARRIER
[157]234  END SUBROUTINE caldyn_BC
[56]235   
[159]236  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
[350]237       f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
[354]238    USE observable_mod
[920]239    USE compute_geopot_mod, ONLY : compute_geopot => compute_geopot_manual
[167]240    USE disvert_mod, ONLY : caldyn_eta, eta_mass
[145]241    USE trace
[129]242    LOGICAL,INTENT(IN)    :: write_out
[126]243    TYPE(t_field),POINTER :: f_phis(:)
[12]244    TYPE(t_field),POINTER :: f_ps(:)
[159]245    TYPE(t_field),POINTER :: f_mass(:)
[126]246    TYPE(t_field),POINTER :: f_theta_rhodz(:)
247    TYPE(t_field),POINTER :: f_u(:)
248    TYPE(t_field),POINTER :: f_q(:)
[350]249    TYPE(t_field),POINTER :: f_geopot(:)
[134]250    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
[360]251    TYPE(t_field) :: f_dps(:)
252    TYPE(t_field) :: f_dmass(:)
253    TYPE(t_field) :: f_dtheta_rhodz(:)
254    TYPE(t_field) :: f_du(:)
[12]255   
[159]256    REAL(rstd),POINTER :: ps(:), dps(:)
[387]257    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:)
[134]258    REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
[159]259    REAL(rstd),POINTER :: qu(:,:)
[162]260    REAL(rstd),POINTER :: qv(:,:)
[151]261
262! temporary shared variable
[529]263    REAL(rstd),POINTER  :: theta(:,:,:) 
[157]264    REAL(rstd),POINTER  :: pk(:,:)
[156]265    REAL(rstd),POINTER  :: geopot(:,:)
[162]266    REAL(rstd),POINTER  :: convm(:,:) 
[151]267    REAL(rstd),POINTER  :: wwuu(:,:)
[157]268       
269    INTEGER :: ind
[151]270    LOGICAL,SAVE :: first=.TRUE.
271!$OMP THREADPRIVATE(first)
[12]272   
[162]273    IF (first) THEN
[151]274      first=.FALSE.
[162]275      IF(caldyn_eta==eta_mass) THEN
276         CALL init_message(f_ps,req_i1,req_ps)
[190]277      ELSE
[162]278         CALL init_message(f_mass,req_i1,req_mass)
279      END IF
[151]280      CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
281      CALL init_message(f_u,req_e1_vect,req_u)
282      CALL init_message(f_qu,req_e1_scal,req_qu)
[361]283      ! Overlapping com/compute (deactivated) : MPI messages need to be sent at first call to caldyn
284      ! This is needed only once : the next ones will be sent by timeloop
[186]285!      IF(caldyn_eta==eta_mass) THEN
286!         CALL send_message(f_ps,req_ps)
287!         CALL wait_message(req_ps) 
288!      ELSE
289!         CALL send_message(f_mass,req_mass)
290!         CALL wait_message(req_mass) 
291!      END IF
292    ENDIF
293   
294    CALL trace_start("caldyn")
295
[361]296    IF(caldyn_eta==eta_mass) THEN
297       CALL send_message(f_ps,req_ps) ! COM00
[529]298       CALL wait_message(req_ps) ! COM00
[361]299    ELSE
300       CALL send_message(f_mass,req_mass) ! COM00
[529]301       CALL wait_message(req_mass) ! COM00
[361]302    END IF
303   
304    CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01
[529]305    CALL wait_message(req_theta_rhodz) ! COM01 Moved from caldyn_pvort
[361]306    CALL send_message(f_u,req_u) ! COM02
[529]307    CALL wait_message(req_u) ! COM02
[126]308
[529]309    IF(.NOT.hydrostatic) THEN
310       STOP 'caldyn_gcm may not be used yet when non-hydrostatic'
311    END IF
312
[126]313    SELECT CASE(caldyn_conserv)
[909]314    CASE(conserv_energy, conserv_gassmann) ! energy-conserving
[128]315       DO ind=1,ndomain
[186]316          IF (.NOT. assigned_domain(ind)) CYCLE
[128]317          CALL swap_dimensions(ind)
318          CALL swap_geometry(ind)
319          ps=f_ps(ind)
[157]320          u=f_u(ind)
321          theta_rhodz = f_theta_rhodz(ind)
[159]322          mass=f_mass(ind)
[157]323          theta = f_theta(ind)
[128]324          qu=f_qu(ind)
[162]325          qv=f_qv(ind)
[529]326          pk = f_pk(ind)
327          geopot = f_geopot(ind) 
328          hflux=f_hflux(ind)
329          convm = f_dmass(ind)
330          dtheta_rhodz=f_dtheta_rhodz(ind)
331          du=f_du(ind)
[387]332          CALL compute_pvort(ps,u,theta_rhodz(:,:,1), mass,theta,qu,qv) ! COM00 COM01 COM02
[529]333!          CALL compute_theta(ps,theta_rhodz, mass,theta)
334!          CALL compute_pvort_only(u,mass,qu,qv)
[128]335
[529]336          CALL compute_geopot(mass,theta, ps,pk,geopot)
337!          du(:,:)=0.
338!          CALL compute_caldyn_fast(0.,u,mass,theta,pk,geopot,du)
339       ENDDO       
[128]340
[529]341       CALL send_message(f_u,req_u) ! COM02
342       CALL wait_message(req_u) ! COM02
343       CALL send_message(f_qu,req_qu) ! COM03
344       CALL wait_message(req_qu) ! COM03
345
[128]346       DO ind=1,ndomain
[186]347          IF (.NOT. assigned_domain(ind)) CYCLE
[128]348          CALL swap_dimensions(ind)
349          CALL swap_geometry(ind)
350          ps=f_ps(ind)
[157]351          u=f_u(ind)
[529]352          theta_rhodz = f_theta_rhodz(ind)
[159]353          mass=f_mass(ind)
[157]354          theta = f_theta(ind)
[128]355          qu=f_qu(ind)
[529]356          qv=f_qv(ind)
[151]357          pk = f_pk(ind)
[156]358          geopot = f_geopot(ind) 
[157]359          hflux=f_hflux(ind)
[162]360          convm = f_dmass(ind)
[157]361          dtheta_rhodz=f_dtheta_rhodz(ind)
362          du=f_du(ind)
[529]363
[387]364          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz(:,:,1),du)
[529]365!          CALL compute_caldyn_slow_hydro(u,mass,hflux,du, .FALSE.) ! FIXME
366!          CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
[162]367          IF(caldyn_eta==eta_mass) THEN
368             wflux=f_wflux(ind)
369             wwuu=f_wwuu(ind)
370             dps=f_dps(ind)
[387]371             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz(:,:,1), du)
[162]372          END IF
[128]373       ENDDO       
[347]374   
[733]375    CASE(conserv_enstrophy) ! enstrophy-conserving
[126]376       DO ind=1,ndomain
[186]377          IF (.NOT. assigned_domain(ind)) CYCLE
[126]378          CALL swap_dimensions(ind)
379          CALL swap_geometry(ind)
[157]380          ps=f_ps(ind)
381          u=f_u(ind)
382          theta_rhodz=f_theta_rhodz(ind)
[159]383          mass=f_mass(ind)
[157]384          theta = f_theta(ind)
385          qu=f_qu(ind)
[182]386          qv=f_qv(ind)
[387]387          CALL compute_pvort(ps,u,theta_rhodz(:,:,1), mass,theta,qu,qv)
[157]388          pk = f_pk(ind)
[156]389          geopot = f_geopot(ind) 
[183]390          CALL compute_geopot(ps,mass,theta, pk,geopot)
[134]391          hflux=f_hflux(ind)
[162]392          convm = f_dmass(ind)
[126]393          dtheta_rhodz=f_dtheta_rhodz(ind)
394          du=f_du(ind)
[387]395          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz(:,:,1),du)
[162]396          IF(caldyn_eta==eta_mass) THEN
397             wflux=f_wflux(ind)
398             wwuu=f_wwuu(ind)
399             dps=f_dps(ind)
400             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
401          END IF
[126]402       ENDDO
403       
404    CASE DEFAULT
[909]405       PRINT *, 'Invalid value of caldyn_conserv : ', caldyn_conserv
[126]406       STOP
407    END SELECT
[12]408
[295]409!$OMP BARRIER
[126]410    !    CALL check_mass_conservation(f_ps,f_dps)
[145]411    CALL trace_end("caldyn")
[186]412!!$OMP BARRIER
[126]413   
414END SUBROUTINE caldyn
[179]415
[126]416!-------------------------------- Diagnostics ----------------------------
417
418  SUBROUTINE check_mass_conservation(f_ps,f_dps)
419    TYPE(t_field),POINTER :: f_ps(:)
420    TYPE(t_field),POINTER :: f_dps(:)
421    REAL(rstd),POINTER :: ps(:)
422    REAL(rstd),POINTER :: dps(:)
423    REAL(rstd) :: mass_tot,dmass_tot
424    INTEGER :: ind,i,j,ij
425   
426    mass_tot=0
427    dmass_tot=0
428   
429    CALL transfert_request(f_dps,req_i1)
430    CALL transfert_request(f_ps,req_i1)
431
432    DO ind=1,ndomain
433      CALL swap_dimensions(ind)
434      CALL swap_geometry(ind)
435
436      ps=f_ps(ind)
437      dps=f_dps(ind)
438
439      DO j=jj_begin,jj_end
440        DO i=ii_begin,ii_end
441          ij=(j-1)*iim+i
442          IF (domain(ind)%own(i,j)) THEN
443            mass_tot=mass_tot+ps(ij)*Ai(ij)/g
444            dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
445          ENDIF
446        ENDDO
447      ENDDO
448   
449    ENDDO
[131]450    IF (is_mpi_root) PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot       
[126]451
452  END SUBROUTINE check_mass_conservation 
[12]453 
454END MODULE caldyn_gcm_mod
Note: See TracBrowser for help on using the repository browser.