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

Last change on this file since 954 was 954, checked in by adurocher, 5 years ago

trunk : Added metric terms to kernels parameters to avoid Host/GPU transferts

Metric terms are now subroutine parameters instead of module variables in kernel subroutines. Dummy arguments for metric terms are now defined as fixed-size arrays, and arrays dimensions are well known when entering an 'acc data' region. Array descriptors are no longer transferred form host to device each time the data region is executed.
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  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    USE abort_mod
18    CHARACTER(len=255) :: def
19    INTEGER            :: ind
20    REAL(rstd),POINTER :: planetvel(:)
21
22    hydrostatic=.TRUE.
23    CALL getin("hydrostatic",hydrostatic)
24
25    IF (.NOT. hydrostatic) THEN
26       CALL abort_acc("hydrostatic /= .TRUE.")
27    END IF
28
29    dysl=.NOT.hydrostatic ! dysl defaults to .FALSE. if hydrostatic, .TRUE. if NH                                             
30    CALL getin("dysl",dysl)
31    dysl_geopot=dysl
32    CALL getin("dysl_geopot",dysl_geopot)
33    dysl_caldyn_fast=dysl
34    CALL getin("dysl_caldyn_fast",dysl_caldyn_fast)
35    dysl_slow_hydro=dysl
36    CALL getin("dysl_slow_hydro",dysl_slow_hydro)
37    dysl_pvort_only=dysl
38    CALL getin("dysl_pvort_only",dysl_pvort_only)
39    dysl_caldyn_coriolis=dysl
40    CALL getin("dysl_caldyn_coriolis",dysl_caldyn_coriolis)
41
42    ! dysl is not supported with openACC
43    IF (dysl_geopot .OR. dysl_caldyn_fast .OR. dysl_slow_hydro .OR. dysl_pvort_only .OR. dysl_caldyn_coriolis) THEN
44       CALL abort_acc("dysl not supported")
45    END IF
46
47    def='energy'
48    CALL getin('caldyn_conserv',def)
49    SELECT CASE(TRIM(def))
50    CASE('energy')
51       caldyn_conserv=energy
52    CASE('enstrophy')
53       caldyn_conserv=enstrophy
54    CASE DEFAULT
55       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', &
56            TRIM(def),'> options are <energy>, <enstrophy>'
57       STOP
58    END SELECT
59    IF (is_master) PRINT *, 'caldyn_conserv=',def
60
61    nqdyn=1 ! default value
62    physics_thermo = thermo_none
63
64    def='theta'
65    CALL getin('thermo',def)
66    SELECT CASE(TRIM(def))
67    CASE('boussinesq')
68       boussinesq=.TRUE.
69       caldyn_thermo=thermo_boussinesq
70       IF(.NOT. hydrostatic) THEN
71          PRINT *, 'thermo=boussinesq and hydrostatic=.FALSE. : Non-hydrostatic boussinesq equations are not supported'
72          STOP
73       END IF
74    CASE('theta')
75       caldyn_thermo=thermo_theta
76       physics_thermo=thermo_dry
77    CASE('entropy')
78       caldyn_thermo=thermo_entropy
79       physics_thermo=thermo_dry
80    CASE('theta_fake_moist')
81       caldyn_thermo=thermo_theta
82       physics_thermo=thermo_fake_moist
83    CASE('entropy_fake_moist')
84       caldyn_thermo=thermo_entropy
85       physics_thermo=thermo_fake_moist
86    CASE('moist')
87       caldyn_thermo=thermo_moist_debug
88       physics_thermo=thermo_moist
89       nqdyn = 2
90    CASE DEFAULT
91       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_thermo : <', &
92            TRIM(def),'> options are <theta>, <entropy>'
93       STOP
94    END SELECT
95
96    IF(is_master) THEN
97       SELECT CASE(caldyn_thermo)
98       CASE(thermo_theta)
99          PRINT *, 'caldyn_thermo = thermo_theta'
100       CASE(thermo_entropy)
101          PRINT *, 'caldyn_thermo = thermo_entropy'
102       CASE(thermo_moist_debug)
103          PRINT *, 'caldyn_thermo = thermo_moist_debug'
104       CASE DEFAULT
105          STOP
106       END SELECT
107
108       SELECT CASE(physics_thermo)
109       CASE(thermo_dry)
110          PRINT *, 'physics_thermo = thermo_dry'
111       CASE(thermo_fake_moist)
112          PRINT *, 'physics_thermo = thermo_fake_moist'
113       CASE(thermo_moist)
114          PRINT *, 'physics_thermo = thermo_moist'
115       END SELECT
116
117       PRINT *, 'nqdyn =', nqdyn
118    END IF
119
120    CALL allocate_caldyn
121
122    DO ind=1,ndomain
123       IF (.NOT. assigned_domain(ind)) CYCLE
124       CALL swap_dimensions(ind)
125       CALL swap_geometry(ind)
126       planetvel=f_planetvel(ind)
127       CALL compute_planetvel(planetvel)
128    END DO
129
130  END SUBROUTINE init_caldyn
131
132  SUBROUTINE allocate_caldyn
133    CALL allocate_field(f_out_u,field_u,type_real,llm) 
134    CALL allocate_field(f_qu,field_u,type_real,llm, ondevice=.TRUE.) 
135    CALL allocate_field(f_qv,field_z,type_real,llm, ondevice=.TRUE.) 
136    CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk', ondevice=.TRUE.)
137    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu', ondevice=.TRUE.)   
138    CALL allocate_field(f_planetvel,  field_u,type_real, name='planetvel') ! planetary velocity at r=a
139    IF(.NOT.hydrostatic) THEN
140       CALL allocate_field(f_Fel,      field_u,type_real,llm+1,name='F_el')
141       CALL allocate_field(f_gradPhi2, field_t,type_real,llm+1,name='gradPhi2')
142       CALL allocate_field(f_wil,      field_t,type_real,llm+1,name='w_il')
143       CALL allocate_field(f_Wetadot,  field_t,type_real,llm,name='W_etadot')
144    END IF
145  END SUBROUTINE allocate_caldyn
146
147  SUBROUTINE caldyn_BC(f_phis, f_geopot, f_wflux)
148    TYPE(t_field),POINTER :: f_phis(:)
149    TYPE(t_field),POINTER :: f_geopot(:)
150    TYPE(t_field),POINTER :: f_wflux(:)
151    REAL(rstd),POINTER  :: phis(:)
152    REAL(rstd),POINTER  :: wflux(:,:)
153    REAL(rstd),POINTER  :: geopot(:,:)
154    REAL(rstd),POINTER  :: wwuu(:,:)
155
156    INTEGER :: ind,ij
157
158    IF (is_omp_first_level) THEN
159       DO ind=1,ndomain
160          IF (.NOT. assigned_domain(ind)) CYCLE
161          CALL swap_dimensions(ind)
162          CALL swap_geometry(ind)
163          geopot=f_geopot(ind)
164          phis=f_phis(ind)
165          wflux=f_wflux(ind)
166          wwuu=f_wwuu(ind)
167         
168          DO ij=ij_begin_ext,ij_end_ext
169              ! lower BCs : geopot=phis, wflux=0, wwuu=0
170              geopot(ij,1) = phis(ij)
171              wflux(ij,1) = 0.
172              wwuu(ij+u_right,1)=0   
173              wwuu(ij+u_lup,1)=0   
174              wwuu(ij+u_ldown,1)=0
175              ! top BCs : wflux=0, wwuu=0
176              wflux(ij,llm+1)  = 0.
177              wwuu(ij+u_right,llm+1)=0   
178              wwuu(ij+u_lup,llm+1)=0   
179              wwuu(ij+u_ldown,llm+1)=0
180          ENDDO
181       END DO
182    ENDIF
183
184    !$OMP BARRIER
185  END SUBROUTINE caldyn_BC
186   
187  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
188       f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
189    USE observable_mod
190    USE disvert_mod, ONLY : caldyn_eta, eta_mass, bp
191    USE trace
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, bp)
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, bp)
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    TYPE(t_field),POINTER :: f_ps(:)
369    TYPE(t_field),POINTER :: f_dps(:)
370    REAL(rstd),POINTER :: ps(:)
371    REAL(rstd),POINTER :: dps(:)
372    REAL(rstd) :: mass_tot,dmass_tot
373    INTEGER :: ind,i,j,ij
374   
375    mass_tot=0
376    dmass_tot=0
377   
378    CALL transfert_request(f_dps,req_i1)
379    CALL transfert_request(f_ps,req_i1)
380
381    DO ind=1,ndomain
382      CALL swap_dimensions(ind)
383      CALL swap_geometry(ind)
384
385      ps=f_ps(ind)
386      dps=f_dps(ind)
387
388      DO j=jj_begin,jj_end
389        DO i=ii_begin,ii_end
390          ij=(j-1)*iim+i
391          IF (domain(ind)%own(i,j)) THEN
392            mass_tot=mass_tot+ps(ij)*Ai(ij)/g
393            dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
394          ENDIF
395        ENDDO
396      ENDDO
397   
398    ENDDO
399    IF (is_mpi_root) PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot       
400
401  END SUBROUTINE check_mass_conservation 
402 
403END MODULE caldyn_gcm_mod
Note: See TracBrowser for help on using the repository browser.