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

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

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