source: codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90 @ 580

Last change on this file since 580 was 580, checked in by dubos, 7 years ago

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