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

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

trunk : GPU implementation with OpenACC ( merge from glcp.idris.fr )
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
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)
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    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.