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

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

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