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

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

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