source: codes/icosagcm/trunk/src/physics/physics.f90 @ 1046

Last change on this file since 1046 was 1046, checked in by ymipsl, 4 years ago

Introduce modification from A. Durocher github to make held&suarez testcase working on GPU

YM & AD
File size: 13.5 KB
Line 
1MODULE physics_mod
2  USE icosa
3  USE field_mod
4  USE physics_interface_mod
5  USE omp_para
6  IMPLICIT NONE
7  PRIVATE
8
9  INTEGER, PARAMETER :: phys_none=0, phys_HS94=1, phys_DCMIP=2, phys_lmdz_generic=3, phys_LB2012=4, phys_external=5,&
10                        phys_DCMIP2016=6
11
12  INTEGER :: phys_type
13  TYPE(t_field),POINTER,SAVE :: f_dulon(:), f_dulat(:)
14  TYPE(t_field),POINTER,SAVE :: f_ulon(:), f_ulat(:)
15  TYPE(t_field),POINTER,SAVE :: f_p(:), f_pk(:)
16  TYPE(t_field),POINTER,SAVE :: f_temp(:)
17  TYPE(t_field),POINTER,SAVE :: f_du_phys(:)
18
19  CHARACTER(LEN=255),SAVE :: physics_type
20!$OMP THREADPRIVATE(physics_type)
21  TYPE(t_message),SAVE :: req_theta0, req_ue0, req_q0
22
23  PUBLIC :: physics, init_physics, zero_du_phys
24
25CONTAINS
26
27  SUBROUTINE init_physics
28    USE mpipara
29    USE abort_mod
30    USE etat0_mod
31    USE physics_dcmip_mod, ONLY : init_physics_dcmip=>init_physics
32    USE physics_dcmip2016_mod, ONLY : init_physics_dcmip2016=>init_physics
33    USE etat0_venus_mod, ONLY : init_phys_venus=>init_physics
34    USE physics_lmdz_generic_mod, ONLY : init_physics_lmdz_generic=>init_physics
35    USE physics_external_mod, ONLY : init_physics_external=>init_physics
36
37    physics_inout%dt_phys = dt*itau_physics
38    physics_type='none'
39    CALL getin("physics",physics_type)
40    SELECT CASE(TRIM(physics_type))
41    CASE ('none')
42
43!$OMP PARALLEL
44       IF(is_mpi_root) PRINT*,"NO PHYSICAL PACKAGE USED"
45       phys_type = phys_none
46!$OMP END PARALLEL
47
48    CASE ('held_suarez')
49
50!$OMP PARALLEL
51       phys_type = phys_HS94
52!$OMP END PARALLEL
53
54    CASE ('Lebonnois2012')
55
56!$OMP PARALLEL
57       phys_type = phys_LB2012
58       CALL init_phys_venus
59!$OMP END PARALLEL
60
61    CASE ('phys_lmdz_generic')
62
63!$OMP PARALLEL
64       CALL init_physics_lmdz_generic
65       phys_type=phys_lmdz_generic
66!$OMP END PARALLEL
67
68    CASE ('phys_external')
69
70       CALL init_physics_external
71!$OMP PARALLEL
72       phys_type=phys_external
73!$OMP END PARALLEL
74
75    CASE ('dcmip')
76
77!$OMP PARALLEL
78       CALL allocate_field(f_dulon,field_t,type_real,llm, name='dulon')
79       CALL allocate_field(f_dulat,field_t,type_real,llm, name='dulat')
80       CALL allocate_field(f_temp,field_t,type_real,llm, name='temp')
81       CALL allocate_field(f_ulon,field_t,type_real,llm, name='ulon')
82       CALL allocate_field(f_ulat,field_t,type_real,llm, name='ulat')
83       CALL allocate_field(f_p,field_t,type_real,llm+1, name='p')
84       CALL allocate_field(f_pk,field_t,type_real,llm, name='pk')
85       CALL init_pack_before ! Compute physics_inout%ngrid and offsets used by pack/unpack
86       CALL init_physics_dcmip
87       CALL init_pack_after ! Defines Ai, lon, lat in physics_inout
88       phys_type = phys_DCMIP
89!$OMP END PARALLEL
90
91    CASE ('dcmip2016')
92
93!$OMP PARALLEL
94       CALL allocate_field(f_dulon,field_t,type_real,llm, name='dulon')
95       CALL allocate_field(f_dulat,field_t,type_real,llm, name='dulat')
96       CALL allocate_field(f_temp,field_t,type_real,llm, name='temp')
97       CALL allocate_field(f_ulon,field_t,type_real,llm, name='ulon')
98       CALL allocate_field(f_ulat,field_t,type_real,llm, name='ulat')
99       CALL allocate_field(f_p,field_t,type_real,llm+1, name='p')
100       CALL allocate_field(f_pk,field_t,type_real,llm, name='pk')
101       CALL init_pack_before ! Compute physics_inout%ngrid and offsets used by pack/unpack
102       CALL init_physics_dcmip2016
103       CALL init_pack_after ! Defines Ai, lon, lat in physics_inout
104       phys_type = phys_DCMIP2016
105!$OMP END PARALLEL
106
107    CASE DEFAULT
108       IF(is_mpi_root) PRINT*, 'init_physics : Bad selector for variable physics <',&
109            TRIM(physics_type), '> options are <none>, <held_suarez>, <Lebonnois2012>, <dcmip>', &
110                                '<phys_lmdz_generic>, <phys_external>'
111       STOP
112    END SELECT
113
114!$OMP PARALLEL
115    CALL allocate_field(f_du_phys,field_u,type_real,llm, name='du_phys', ondevice=.TRUE.)
116
117    IF(is_mpi_root) PRINT *, 'phys_type = ',phys_type
118   
119   
120!$OMP END PARALLEL
121    IF (phys_type /= phys_none .AND. phys_type /= phys_HS94) THEN
122       CALL abort_acc("only phys_type == 'phys_none' or 'phys_HS94' supported")
123    END IF
124  END SUBROUTINE init_physics
125
126  SUBROUTINE zero_du_phys()
127    USE abort_mod
128    REAL(rstd), DIMENSION(:,:), POINTER :: du
129    INTEGER :: ind
130    DO ind=1,ndomain
131       IF (.NOT. assigned_domain(ind)) CYCLE
132       CALL swap_dimensions(ind)
133       CALL swap_geometry(ind)
134       du=f_du_phys(ind)
135
136       !$acc kernels default(present) async
137       du(:,ll_begin:ll_end) = 0.
138       !$acc end kernels
139         
140    END DO
141  END SUBROUTINE zero_du_phys
142
143  SUBROUTINE add_du_phys(coef, f_u)
144    REAL(rstd), INTENT(IN) :: coef  ! -1 before physics, +1 after physics
145    TYPE(t_field),POINTER :: f_u(:) ! velocity field before/after call to physics
146    REAL(rstd), DIMENSION(:,:), POINTER :: u, du
147    INTEGER :: ind, ij
148
149    DO ind=1,ndomain
150       IF (.NOT. assigned_domain(ind)) CYCLE
151       CALL swap_dimensions(ind)
152       CALL swap_geometry(ind)
153       du=f_du_phys(ind)
154       u=f_u(ind)
155
156       !$acc parallel loop default(present) async
157       DO ij = ij_begin, ij_end
158          du(ij,ll_begin:ll_end) = du(ij,ll_begin:ll_end) + coef*u(ij,ll_begin:ll_end)
159       END DO
160
161    END DO
162  END SUBROUTINE add_du_phys
163
164  SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
165    USE physics_lmdz_generic_mod, ONLY : physics_lmdz_generic => physics
166    USE physics_external_mod, ONLY : physics_external => physics
167    USE physics_dcmip_mod, ONLY : write_physics_dcmip => write_physics
168    USE physics_dcmip2016_mod, ONLY : write_physics_dcmip2016 => write_physics
169    USE etat0_heldsz_mod
170    USE etat0_venus_mod, ONLY : phys_venus => physics
171    USE abort_mod
172    INTEGER, INTENT(IN)   :: it
173    TYPE(t_field),POINTER :: f_phis(:)
174    TYPE(t_field),POINTER :: f_ps(:)
175    TYPE(t_field),POINTER :: f_theta_rhodz(:)
176    TYPE(t_field),POINTER :: f_ue(:)
177    TYPE(t_field),POINTER :: f_wflux(:)
178    TYPE(t_field),POINTER :: f_q(:)
179
180    LOGICAL,SAVE :: first=.TRUE.
181    !$OMP THREADPRIVATE(first)
182    IF (first) THEN
183      CALL init_message(f_theta_rhodz, req_i0, req_theta0)
184      CALL init_message(f_ue, req_e0_vect, req_ue0)
185      CALL init_message(f_q, req_i0, req_q0)
186      first=.FALSE.
187    ENDIF
188   
189    IF (phys_type==phys_external) THEN
190   
191      CALL physics_external(it ,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
192
193    ELSE
194   
195      IF(MOD(it,itau_physics)==0 .AND. phys_type/=phys_none) THEN
196
197         ! as a result of the the two calls to add_du_phys,
198         ! du_phys increases by u(after physics) - u (before physics)
199         CALL add_du_phys(-1., f_ue)
200
201         SELECT CASE(phys_type)
202         CASE(phys_HS94)
203            CALL held_suarez(f_ps,f_theta_rhodz,f_ue) 
204         CASE (phys_lmdz_generic)
205           CALL physics_lmdz_generic(it ,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
206         CASE(phys_LB2012)
207            CALL phys_venus(f_ps,f_theta_rhodz,f_ue) 
208         CASE DEFAULT
209            CALL physics_column(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
210         END SELECT
211
212         CALL send_message(f_theta_rhodz, req_theta0)
213         CALL send_message(f_ue, req_ue0)
214         CALL send_message(f_q, req_q0)
215         CALL wait_message(req_theta0)
216         CALL wait_message(req_ue0)
217         CALL wait_message(req_q0)
218       
219         CALL add_du_phys(1., f_ue)
220      END IF
221
222      IF (mod(it,itau_out)==0 ) THEN
223         CALL write_physics_tendencies
224         CALL zero_du_phys
225         SELECT CASE(phys_type)
226         CASE (phys_DCMIP)
227            CALL write_physics_dcmip
228         CASE (phys_DCMIP2016)
229            CALL write_physics_dcmip2016
230         END SELECT
231      END IF
232    ENDIF
233   
234  END SUBROUTINE physics
235
236  SUBROUTINE write_physics_tendencies
237    USE observable_mod, ONLY : f_buf_ulon, f_buf_ulat
238    USE wind_mod
239    USE output_field_mod
240
241    CALL transfert_request(f_du_phys,req_e1_vect)   
242    CALL update_host_field(f_du_phys)
243
244    CALL un2ulonlat(f_du_phys, f_buf_ulon, f_buf_ulat, (1./(dt*itau_out)))
245    CALL output_field("dulon_phys",f_buf_ulon)
246    CALL output_field("dulat_phys",f_buf_ulat)
247  END SUBROUTINE write_physics_tendencies
248   
249  SUBROUTINE physics_column(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
250    USE physics_dcmip_mod, ONLY : full_physics_dcmip => full_physics
251    USE physics_dcmip2016_mod, ONLY : full_physics_dcmip2016 => full_physics
252    USE theta2theta_rhodz_mod
253    USE mpipara
254    USE checksum_mod
255    TYPE(t_field),POINTER :: f_phis(:)
256    TYPE(t_field),POINTER :: f_ps(:)
257    TYPE(t_field),POINTER :: f_theta_rhodz(:)
258    TYPE(t_field),POINTER :: f_ue(:)
259    TYPE(t_field),POINTER :: f_q(:)
260    REAL(rstd),POINTER :: phis(:)
261    REAL(rstd),POINTER :: ps(:)
262    REAL(rstd),POINTER :: temp(:,:)
263    REAL(rstd),POINTER :: ue(:,:)
264    REAL(rstd),POINTER :: dulon(:,:)
265    REAL(rstd),POINTER :: dulat(:,:)
266    REAL(rstd),POINTER :: q(:,:,:)
267    REAL(rstd),POINTER :: p(:,:)
268    REAL(rstd),POINTER :: pk(:,:)
269    REAL(rstd),POINTER :: ulon(:,:)
270    REAL(rstd),POINTER :: ulat(:,:)
271    INTEGER :: it, ind
272
273    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_temp)
274   
275    DO ind=1,ndomain
276       IF (.NOT. assigned_domain(ind)) CYCLE
277       CALL swap_dimensions(ind)
278       CALL swap_geometry(ind)
279       phis=f_phis(ind)
280       ps=f_ps(ind)
281       temp=f_temp(ind)
282       ue=f_ue(ind)
283       q=f_q(ind)
284       p=f_p(ind)
285       pk=f_pk(ind)
286       ulon=f_ulon(ind)
287       ulat=f_ulat(ind)
288       CALL pack_physics(pack_info(ind), phis, ps, temp, ue, q, p, pk, ulon, ulat)
289    END DO
290
291    SELECT CASE(phys_type)
292    CASE (phys_DCMIP)
293       IF (is_omp_level_master) CALL full_physics_dcmip
294    CASE (phys_DCMIP2016)
295       IF (is_omp_level_master) CALL full_physics_dcmip2016
296    CASE DEFAULT
297       IF(is_master) PRINT *,'Internal error : illegal value of phys_type', phys_type
298       STOP
299    END SELECT
300
301    DO ind=1,ndomain
302       IF (.NOT. assigned_domain(ind)) CYCLE
303       CALL swap_dimensions(ind)
304       CALL swap_geometry(ind)
305       ps=f_ps(ind)
306       temp=f_temp(ind)
307       q=f_q(ind)
308       dulon=f_dulon(ind)
309       dulat=f_dulat(ind)
310       CALL unpack_physics(pack_info(ind), ps, temp, q, dulon, dulat)
311    END DO
312   
313    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
314
315    ! Transfer dulon, dulat
316    CALL transfert_request(f_dulon,req_i0)
317    CALL transfert_request(f_dulat,req_i0)
318
319    DO ind=1,ndomain
320       IF (.NOT. assigned_domain(ind)) CYCLE
321       CALL swap_dimensions(ind)
322       CALL swap_geometry(ind)
323       ue=f_ue(ind)
324       dulon=f_dulon(ind)
325       dulat=f_dulat(ind)
326       CALL compute_update_velocity(dulon, dulat, ue)
327    END DO
328
329  END SUBROUTINE physics_column
330
331  SUBROUTINE pack_physics(info, phis, ps, temp, ue, q, p, pk, ulon, ulat )
332    USE wind_mod
333    USE pression_mod
334    USE theta2theta_rhodz_mod
335    USE exner_mod
336    TYPE(t_pack_info) :: info
337    REAL(rstd) :: phis(iim*jjm)
338    REAL(rstd) :: ps(iim*jjm)
339    REAL(rstd) :: temp(iim*jjm,llm)
340    REAL(rstd) :: pks(iim*jjm)
341    REAL(rstd) :: pk(iim*jjm,llm)
342    REAL(rstd) :: ue(3*iim*jjm,llm)
343    REAL(rstd) :: q(iim*jjm,llm,nqtot)
344
345    REAL(rstd) :: p(iim*jjm,llm+1)
346    REAL(rstd) :: uc(iim*jjm,llm,3)
347    REAL(rstd) :: ulon(iim*jjm,llm)
348    REAL(rstd) :: ulat(iim*jjm,llm)
349
350!$OMP BARRIER
351    CALL compute_pression(ps,p,0)
352!$OMP BARRIER
353    CALL compute_exner(ps,p,pks,pk,0) 
354!$OMP BARRIER
355    CALL compute_wind_centered(ue,uc)
356    CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
357!$OMP BARRIER
358    IF (is_omp_level_master) THEN
359      CALL pack_domain(info, phis, physics_inout%phis)
360      CALL pack_domain(info, p, physics_inout%p)
361      CALL pack_domain(info, pk, physics_inout%pk)
362      CALL pack_domain(info, Temp, physics_inout%Temp)
363      CALL pack_domain(info, ulon, physics_inout%ulon)
364      CALL pack_domain(info, ulat, physics_inout%ulat)
365      CALL pack_domain(info, q, physics_inout%q)
366    ENDIF
367!$OMP BARRIER
368  END SUBROUTINE pack_physics
369
370  SUBROUTINE unpack_physics(info, ps,temp, q, dulon, dulat)
371    USE theta2theta_rhodz_mod
372    TYPE(t_pack_info) :: info
373    REAL(rstd) :: ps(iim*jjm)
374    REAL(rstd) :: temp(iim*jjm,llm)
375    REAL(rstd) :: q(iim*jjm,llm,nqtot)
376    REAL(rstd) :: dulon(iim*jjm,llm)
377    REAL(rstd) :: dulat(iim*jjm,llm)
378
379    REAL(rstd) :: dq(iim*jjm,llm,nqtot)
380    REAL(rstd) :: dTemp(iim*jjm,llm)
381
382!$OMP BARRIER
383    IF (is_omp_level_master) THEN
384      CALL unpack_domain(info, dulon, physics_inout%dulon)
385      CALL unpack_domain(info, dulat, physics_inout%dulat)
386      CALL unpack_domain(info, dq, physics_inout%dq)
387      CALL unpack_domain(info, Temp, physics_inout%Temp)
388      CALL unpack_domain(info, dTemp, physics_inout%dTemp)
389      q = q + physics_inout%dt_phys * dq
390      Temp = Temp + physics_inout%dt_phys * dTemp
391    ENDIF
392!$OMP BARRIER
393
394!    CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
395  END SUBROUTINE unpack_physics
396
397  SUBROUTINE compute_update_velocity(dulon, dulat, ue)
398    USE wind_mod
399    REAL(rstd) :: dulon(iim*jjm,llm)
400    REAL(rstd) :: dulat(iim*jjm,llm)
401    REAL(rstd) :: ue(3*iim*jjm,llm)
402    REAL(rstd) :: duc(iim*jjm,llm,3)
403    REAL(rstd) :: dt2, due
404    INTEGER :: i,j,ij,l
405    ! Reconstruct wind tendencies at edges and add to normal wind
406    CALL compute_wind_centered_from_lonlat_compound(dulon,dulat,duc)
407    dt2=.5*physics_inout%dt_phys
408    DO l=ll_begin,ll_end
409      DO j=jj_begin,jj_end
410        DO i=ii_begin,ii_end
411          ij=(j-1)*iim+i
412          due = sum( (duc(ij,l,:) + duc(ij+t_right,l,:))*ep_e(ij+u_right,:) )
413          ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due
414
415          due = sum( (duc(ij,l,:) + duc(ij+t_lup,l,:))*ep_e(ij+u_lup,:) )
416          ue(ij+u_lup,l)=ue(ij+u_lup,l) + dt2*due
417
418          due = sum( (duc(ij,l,:) + duc(ij+t_ldown,l,:))*ep_e(ij+u_ldown,:) )
419          ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + dt2*due
420        ENDDO
421      ENDDO
422    ENDDO
423  END SUBROUTINE compute_update_velocity
424
425END MODULE physics_mod
Note: See TracBrowser for help on using the repository browser.