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

Last change on this file since 381 was 381, checked in by ymipsl, 8 years ago

Add DCMIP2016 physics

(first guess)

YM
File size: 10.2 KB
Line 
1MODULE physics_mod
2
3  USE field_mod
4
5  PRIVATE
6
7  INTEGER, PARAMETER :: phys_none=0, phys_HS94=1, phys_DCMIP=2, phys_lmdz_generic=3, phys_LB2012=4, phys_external=5, phys_DCMIP2016=6
8
9  INTEGER :: phys_type
10  TYPE(t_field),POINTER :: f_extra_physics_2D(:), f_extra_physics_3D(:)
11  TYPE(t_field),POINTER :: f_dulon(:), f_dulat(:)
12  TYPE(t_field),POINTER :: f_temp(:)
13
14  CHARACTER(LEN=255) :: physics_type
15!$OMP THREADPRIVATE(physics_type)
16
17  PUBLIC :: physics, init_physics
18
19CONTAINS
20
21  SUBROUTINE init_physics
22    USE mpipara
23    USE etat0_mod
24    USE icosa
25    USE physics_interface_mod
26    USE physics_dcmip_mod, ONLY : init_physics_dcmip=>init_physics
27    USE physics_dcmip2016_mod, ONLY : init_physics_dcmip2016=>init_physics
28    USE etat0_venus_mod, ONLY : init_phys_venus=>init_physics
29    USE physics_lmdz_generic_mod, ONLY : init_physics_lmdz_generic=>init_physics
30    USE physics_external_mod, ONLY : init_physics_external=>init_physics
31    IMPLICIT NONE
32
33    physics_inout%dt_phys = dt*itau_physics
34    physics_type='none'
35    CALL getin("physics",physics_type)
36    SELECT CASE(TRIM(physics_type))
37    CASE ('none')
38       IF(is_mpi_root) PRINT*,"NO PHYSICAL PACKAGE USED"
39       phys_type = phys_none
40    CASE ('held_suarez')
41       phys_type = phys_HS94
42    CASE ('Lebonnois2012')
43       phys_type = phys_LB2012
44       CALL init_phys_venus
45
46    CASE ('phys_lmdz_generic')
47       CALL init_physics_lmdz_generic
48       phys_type=phys_lmdz_generic
49    CASE ('phys_external')
50       CALL init_physics_external
51       phys_type=phys_external
52    CASE ('dcmip')
53       CALL allocate_field(f_dulon,field_t,type_real,llm, name='dulon')
54       CALL allocate_field(f_dulat,field_t,type_real,llm, name='dulat')
55       CALL allocate_field(f_temp,field_t,type_real,llm, name='temp')
56       CALL init_pack_before ! Compute physics_inout%ngrid and offsets used by pack/unpack
57       CALL init_physics_dcmip
58       CALL init_pack_after ! Defines Ai, lon, lat in physics_inout
59       phys_type = phys_DCMIP
60    CASE ('dcmip2016')
61       CALL allocate_field(f_dulon,field_t,type_real,llm, name='dulon')
62       CALL allocate_field(f_dulat,field_t,type_real,llm, name='dulat')
63       CALL allocate_field(f_temp,field_t,type_real,llm, name='temp')
64       CALL allocate_field(f_temp,field_t,type_real,llm, name='temp')
65       CALL init_pack_before ! Compute physics_inout%ngrid and offsets used by pack/unpack
66       CALL init_physics_dcmip2016
67       CALL init_pack_after ! Defines Ai, lon, lat in physics_inout
68       phys_type = phys_DCMIP2016
69    CASE DEFAULT
70       IF(is_mpi_root) PRINT*, 'init_physics : Bad selector for variable physics <',&
71            TRIM(physics_type), '> options are <none>, <held_suarez>, <Lebonnois2012>, <dcmip>', &
72                                '<phys_lmdz_generic>, <phys_external>'
73       STOP
74    END SELECT
75
76    IF(is_mpi_root) PRINT *, 'phys_type = ',phys_type
77  END SUBROUTINE init_physics
78
79  SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
80    USE icosa
81    USE physics_interface_mod
82    USE physics_lmdz_generic_mod, ONLY : physics_lmdz_generic => physics
83    USE physics_external_mod, ONLY : physics_external => physics
84    USE physics_dcmip_mod, ONLY : write_physics_dcmip => write_physics
85    USE physics_dcmip2016_mod, ONLY : write_physics_dcmip2016 => write_physics
86    USE etat0_heldsz_mod
87    USE etat0_venus_mod, ONLY : phys_venus => physics
88    IMPLICIT NONE
89    INTEGER, INTENT(IN)   :: it
90    TYPE(t_field),POINTER :: f_phis(:)
91    TYPE(t_field),POINTER :: f_ps(:)
92    TYPE(t_field),POINTER :: f_theta_rhodz(:)
93    TYPE(t_field),POINTER :: f_ue(:)
94    TYPE(t_field),POINTER :: f_wflux(:)
95    TYPE(t_field),POINTER :: f_q(:)
96    REAL(rstd),POINTER :: phis(:)
97    REAL(rstd),POINTER :: ps(:)
98    REAL(rstd),POINTER :: theta_rhodz(:,:)
99    REAL(rstd),POINTER :: ue(:,:)
100    REAL(rstd),POINTER :: q(:,:,:)
101
102    LOGICAL:: firstcall,lastcall
103    INTEGER :: ind
104    TYPE(t_physics_inout) :: args
105
106    IF(MOD(it,itau_physics)==0) THEN
107   
108       SELECT CASE(phys_type)
109       CASE (phys_none)
110          ! No physics, do nothing
111       CASE(phys_HS94)
112          CALL held_suarez(f_ps,f_theta_rhodz,f_ue) 
113       CASE (phys_lmdz_generic)
114         CALL physics_lmdz_generic(it ,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
115       CASE (phys_external)
116         CALL physics_external(it ,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
117       CASE(phys_LB2012)
118          CALL phys_venus(f_ps,f_theta_rhodz,f_ue) 
119       CASE DEFAULT
120          CALL physics_column(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
121       END SELECT
122
123       CALL transfert_request(f_theta_rhodz,req_i0)
124       CALL transfert_request(f_ue,req_e0_vect)
125       CALL transfert_request(f_q,req_i0)
126    END IF
127
128    IF (mod(it,itau_out)==0 ) THEN
129       SELECT CASE(phys_type)
130       CASE (phys_DCMIP)
131          CALL write_physics_dcmip
132       CASE (phys_DCMIP2016)
133          CALL write_physics_dcmip2016
134       END SELECT
135    END IF
136   
137  END SUBROUTINE physics
138
139  SUBROUTINE physics_column(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
140    USE icosa
141    USE physics_interface_mod
142    USE physics_dcmip_mod, ONLY : full_physics_dcmip => full_physics
143    USE physics_dcmip2016_mod, ONLY : full_physics_dcmip2016 => full_physics
144    USE theta2theta_rhodz_mod
145    USE mpipara
146    IMPLICIT NONE
147    TYPE(t_field),POINTER :: f_phis(:)
148    TYPE(t_field),POINTER :: f_ps(:)
149    TYPE(t_field),POINTER :: f_theta_rhodz(:)
150    TYPE(t_field),POINTER :: f_ue(:)
151    TYPE(t_field),POINTER :: f_q(:)
152    REAL(rstd),POINTER :: phis(:)
153    REAL(rstd),POINTER :: ps(:)
154    REAL(rstd),POINTER :: temp(:,:)
155    REAL(rstd),POINTER :: theta_rhodz(:,:)
156    REAL(rstd),POINTER :: ue(:,:)
157    REAL(rstd),POINTER :: dulon(:,:)
158    REAL(rstd),POINTER :: dulat(:,:)
159    REAL(rstd),POINTER :: q(:,:,:)
160    INTEGER :: it, ind
161
162    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_temp)
163   
164    DO ind=1,ndomain
165       IF (.NOT. assigned_domain(ind)) CYCLE
166       CALL swap_dimensions(ind)
167       CALL swap_geometry(ind)
168       phis=f_phis(ind)
169       ps=f_ps(ind)
170       temp=f_temp(ind)
171       ue=f_ue(ind)
172       q=f_q(ind)
173       CALL pack_physics(pack_info(ind), phis, ps, temp, ue, q)
174    END DO
175
176    SELECT CASE(phys_type)
177    CASE (phys_DCMIP)
178       CALL full_physics_dcmip
179    CASE (phys_DCMIP2016)
180       CALL full_physics_dcmip2016
181    CASE DEFAULT
182       IF(is_mpi_master) PRINT *,'Internal error : illegal value of phys_type', phys_type
183       STOP
184    END SELECT
185
186    DO ind=1,ndomain
187       IF (.NOT. assigned_domain(ind)) CYCLE
188       CALL swap_dimensions(ind)
189       CALL swap_geometry(ind)
190       ps=f_ps(ind)
191       temp=f_temp(ind)
192       q=f_q(ind)
193       dulon=f_dulon(ind)
194       dulat=f_dulat(ind)
195       CALL unpack_physics(pack_info(ind), ps, temp, q, dulon, dulat)
196    END DO
197    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
198
199    ! Transfer dulon, dulat
200    CALL transfert_request(f_dulon,req_i0)
201    CALL transfert_request(f_dulat,req_i0)
202
203    DO ind=1,ndomain
204       IF (.NOT. assigned_domain(ind)) CYCLE
205       CALL swap_dimensions(ind)
206       CALL swap_geometry(ind)
207       ue=f_ue(ind)
208       dulon=f_dulon(ind)
209       dulat=f_dulat(ind)
210       CALL compute_update_velocity(dulon, dulat, ue)
211    END DO
212
213  END SUBROUTINE physics_column
214
215  SUBROUTINE pack_physics(info, phis, ps, temp, ue, q )
216    USE icosa
217    USE wind_mod
218    USE pression_mod
219    USE theta2theta_rhodz_mod
220    USE physics_interface_mod
221    USE exner_mod
222    IMPLICIT NONE
223    TYPE(t_pack_info) :: info
224    REAL(rstd) :: phis(iim*jjm)
225    REAL(rstd) :: ps(iim*jjm)
226    REAL(rstd) :: temp(iim*jjm,llm)
227    REAL(rstd) :: pks(iim*jjm)
228    REAL(rstd) :: pk(iim*jjm,llm)
229    REAL(rstd) :: ue(3*iim*jjm,llm)
230    REAL(rstd) :: q(iim*jjm,llm,nqtot)
231
232    REAL(rstd) :: p(iim*jjm,llm+1)
233    REAL(rstd) :: uc(iim*jjm,3,llm)
234    REAL(rstd) :: ulon(iim*jjm,llm)
235    REAL(rstd) :: ulat(iim*jjm,llm)
236
237!$OMP BARRIER
238    CALL compute_pression(ps,p,0)
239!$OMP BARRIER
240    CALL compute_exner(ps,p,pks,pk,0) 
241!$OMP BARRIER
242    CALL compute_wind_centered(ue,uc)
243    CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
244
245    CALL pack_domain(info, phis, physics_inout%phis)
246    CALL pack_domain(info, p, physics_inout%p)
247    CALL pack_domain(info, pk, physics_inout%pk)
248    CALL pack_domain(info, Temp, physics_inout%Temp)
249    CALL pack_domain(info, ulon, physics_inout%ulon)
250    CALL pack_domain(info, ulat, physics_inout%ulat)
251    CALL pack_domain(info, q, physics_inout%q)
252  END SUBROUTINE pack_physics
253
254  SUBROUTINE unpack_physics(info, ps,temp, q, dulon, dulat)
255    USE icosa
256    USE physics_interface_mod
257    USE theta2theta_rhodz_mod
258    IMPLICIT NONE
259    TYPE(t_pack_info) :: info
260    REAL(rstd) :: ps(iim*jjm)
261    REAL(rstd) :: temp(iim*jjm,llm)
262    REAL(rstd) :: q(iim*jjm,llm,nqtot)
263    REAL(rstd) :: dulon(iim*jjm,llm)
264    REAL(rstd) :: dulat(iim*jjm,llm)
265
266    REAL(rstd) :: dq(iim*jjm,llm,nqtot)
267    REAL(rstd) :: dTemp(iim*jjm,llm)
268    CALL unpack_domain(info, dulon, physics_inout%dulon)
269    CALL unpack_domain(info, dulat, physics_inout%dulat)
270    CALL unpack_domain(info, dq, physics_inout%dq)
271    CALL unpack_domain(info, Temp, physics_inout%Temp)
272    CALL unpack_domain(info, dTemp, physics_inout%dTemp)
273    q = q + physics_inout%dt_phys * dq
274    Temp = Temp + physics_inout%dt_phys * dTemp
275!    CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
276  END SUBROUTINE unpack_physics
277
278  SUBROUTINE compute_update_velocity(dulon, dulat, ue)
279    USE icosa
280    USE physics_interface_mod
281    USE wind_mod
282    IMPLICIT NONE
283    REAL(rstd) :: dulon(iim*jjm,llm)
284    REAL(rstd) :: dulat(iim*jjm,llm)
285    REAL(rstd) :: ue(3*iim*jjm,llm)
286    REAL(rstd) :: duc(iim*jjm,3,llm)
287    REAL(rstd) :: dt2, due
288    INTEGER :: i,j,ij,l
289    ! Reconstruct wind tendencies at edges and add to normal wind
290    CALL compute_wind_centered_from_lonlat_compound(dulon,dulat,duc)
291    dt2=.5*physics_inout%dt_phys
292    DO l=1,llm
293      DO j=jj_begin,jj_end
294        DO i=ii_begin,ii_end
295          ij=(j-1)*iim+i
296          due = sum( (duc(ij,:,l) + duc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
297          ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due
298
299          due = sum( (duc(ij,:,l) + duc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )
300          ue(ij+u_lup,l)=ue(ij+u_lup,l) + dt2*due
301
302          due = sum( (duc(ij,:,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )
303          ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + dt2*due
304        ENDDO
305      ENDDO
306    ENDDO
307  END SUBROUTINE compute_update_velocity
308
309END MODULE physics_mod
Note: See TracBrowser for help on using the repository browser.