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

Last change on this file since 714 was 714, checked in by dubos, 6 years ago

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