source: codes/icosagcm/devel/src/physics/physics.f90

Last change on this file was 1029, checked in by dubos, 4 years ago

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