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

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

bug fix when porting to gortran (mostly line size >132)

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