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

Last change on this file since 584 was 584, checked in by dubos, 7 years ago

devel : switched indices of 3D vector fields
File size: 11.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,&
8                        phys_DCMIP2016=6
9
10  INTEGER :: phys_type
11  TYPE(t_field),POINTER,SAVE :: f_extra_physics_2D(:), f_extra_physics_3D(:)
12  TYPE(t_field),POINTER,SAVE :: f_dulon(:), f_dulat(:)
13  TYPE(t_field),POINTER,SAVE :: f_ulon(:), f_ulat(:)
14  TYPE(t_field),POINTER,SAVE :: f_p(:), f_pk(:)
15  TYPE(t_field),POINTER,SAVE :: f_temp(:)
16
17  CHARACTER(LEN=255),SAVE :: physics_type
18!$OMP THREADPRIVATE(physics_type)
19
20  PUBLIC :: physics, init_physics
21
22CONTAINS
23
24  SUBROUTINE init_physics
25    USE mpipara
26    USE etat0_mod
27    USE icosa
28    USE physics_interface_mod
29    USE physics_dcmip_mod, ONLY : init_physics_dcmip=>init_physics
30    USE physics_dcmip2016_mod, ONLY : init_physics_dcmip2016=>init_physics
31    USE etat0_venus_mod, ONLY : init_phys_venus=>init_physics
32    USE physics_lmdz_generic_mod, ONLY : init_physics_lmdz_generic=>init_physics
33    USE physics_external_mod, ONLY : init_physics_external=>init_physics
34    IMPLICIT NONE
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    IF(is_mpi_root) PRINT *, 'phys_type = ',phys_type
87  END SUBROUTINE init_physics
88
89  SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
90    USE icosa
91    USE physics_interface_mod
92    USE physics_lmdz_generic_mod, ONLY : physics_lmdz_generic => physics
93    USE physics_external_mod, ONLY : physics_external => physics
94    USE physics_dcmip_mod, ONLY : write_physics_dcmip => write_physics
95    USE physics_dcmip2016_mod, ONLY : write_physics_dcmip2016 => write_physics
96    USE etat0_heldsz_mod
97    USE etat0_venus_mod, ONLY : phys_venus => physics
98    IMPLICIT NONE
99    INTEGER, INTENT(IN)   :: it
100    TYPE(t_field),POINTER :: f_phis(:)
101    TYPE(t_field),POINTER :: f_ps(:)
102    TYPE(t_field),POINTER :: f_theta_rhodz(:)
103    TYPE(t_field),POINTER :: f_ue(:)
104    TYPE(t_field),POINTER :: f_wflux(:)
105    TYPE(t_field),POINTER :: f_q(:)
106
107    LOGICAL:: firstcall,lastcall
108    INTEGER :: ind
109    TYPE(t_physics_inout) :: args
110
111    IF(MOD(it,itau_physics)==0) THEN
112   
113       SELECT CASE(phys_type)
114       CASE (phys_none)
115          ! No physics, do nothing
116       CASE(phys_HS94)
117          CALL held_suarez(f_ps,f_theta_rhodz,f_ue) 
118       CASE (phys_lmdz_generic)
119         CALL physics_lmdz_generic(it ,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
120       CASE (phys_external)
121         CALL physics_external(it ,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
122       CASE(phys_LB2012)
123          CALL phys_venus(f_ps,f_theta_rhodz,f_ue) 
124       CASE DEFAULT
125          CALL physics_column(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
126       END SELECT
127
128       CALL transfert_request(f_theta_rhodz,req_i0)
129       CALL transfert_request(f_ue,req_e0_vect)
130       CALL transfert_request(f_q,req_i0)
131    END IF
132
133    IF (mod(it,itau_out)==0 ) THEN
134       SELECT CASE(phys_type)
135       CASE (phys_DCMIP)
136          CALL write_physics_dcmip
137       CASE (phys_DCMIP2016)
138          CALL write_physics_dcmip2016
139       END SELECT
140    END IF
141   
142  END SUBROUTINE physics
143
144  SUBROUTINE physics_column(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
145    USE icosa
146    USE physics_interface_mod
147    USE physics_dcmip_mod, ONLY : full_physics_dcmip => full_physics
148    USE physics_dcmip2016_mod, ONLY : full_physics_dcmip2016 => full_physics
149    USE theta2theta_rhodz_mod
150    USE mpipara
151    USE omp_para
152    USE checksum_mod
153    IMPLICIT NONE
154    TYPE(t_field),POINTER :: f_phis(:)
155    TYPE(t_field),POINTER :: f_ps(:)
156    TYPE(t_field),POINTER :: f_theta_rhodz(:)
157    TYPE(t_field),POINTER :: f_ue(:)
158    TYPE(t_field),POINTER :: f_q(:)
159    REAL(rstd),POINTER :: phis(:)
160    REAL(rstd),POINTER :: ps(:)
161    REAL(rstd),POINTER :: temp(:,:)
162    REAL(rstd),POINTER :: ue(:,:)
163    REAL(rstd),POINTER :: dulon(:,:)
164    REAL(rstd),POINTER :: dulat(:,:)
165    REAL(rstd),POINTER :: q(:,:,:)
166    REAL(rstd),POINTER :: p(:,:)
167    REAL(rstd),POINTER :: pk(:,:)
168    REAL(rstd),POINTER :: ulon(:,:)
169    REAL(rstd),POINTER :: ulat(:,:)
170    INTEGER :: it, ind
171
172    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_temp)
173   
174    DO ind=1,ndomain
175       IF (.NOT. assigned_domain(ind)) CYCLE
176       CALL swap_dimensions(ind)
177       CALL swap_geometry(ind)
178       phis=f_phis(ind)
179       ps=f_ps(ind)
180       temp=f_temp(ind)
181       ue=f_ue(ind)
182       q=f_q(ind)
183       p=f_p(ind)
184       pk=f_pk(ind)
185       ulon=f_ulon(ind)
186       ulat=f_ulat(ind)
187       CALL pack_physics(pack_info(ind), phis, ps, temp, ue, q, p, pk, ulon, ulat)
188    END DO
189
190    SELECT CASE(phys_type)
191    CASE (phys_DCMIP)
192       IF (is_omp_level_master) CALL full_physics_dcmip
193    CASE (phys_DCMIP2016)
194       IF (is_omp_level_master) CALL full_physics_dcmip2016
195    CASE DEFAULT
196       IF(is_master) PRINT *,'Internal error : illegal value of phys_type', phys_type
197       STOP
198    END SELECT
199
200    DO ind=1,ndomain
201       IF (.NOT. assigned_domain(ind)) CYCLE
202       CALL swap_dimensions(ind)
203       CALL swap_geometry(ind)
204       ps=f_ps(ind)
205       temp=f_temp(ind)
206       q=f_q(ind)
207       dulon=f_dulon(ind)
208       dulat=f_dulat(ind)
209       CALL unpack_physics(pack_info(ind), ps, temp, q, dulon, dulat)
210    END DO
211   
212    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
213
214    ! Transfer dulon, dulat
215    CALL transfert_request(f_dulon,req_i0)
216    CALL transfert_request(f_dulat,req_i0)
217
218    DO ind=1,ndomain
219       IF (.NOT. assigned_domain(ind)) CYCLE
220       CALL swap_dimensions(ind)
221       CALL swap_geometry(ind)
222       ue=f_ue(ind)
223       dulon=f_dulon(ind)
224       dulat=f_dulat(ind)
225       CALL compute_update_velocity(dulon, dulat, ue)
226    END DO
227
228  END SUBROUTINE physics_column
229
230  SUBROUTINE pack_physics(info, phis, ps, temp, ue, q, p, pk, ulon, ulat )
231    USE icosa
232    USE wind_mod
233    USE pression_mod
234    USE theta2theta_rhodz_mod
235    USE physics_interface_mod
236    USE exner_mod
237    USE omp_para
238    IMPLICIT NONE
239    TYPE(t_pack_info) :: info
240    REAL(rstd) :: phis(iim*jjm)
241    REAL(rstd) :: ps(iim*jjm)
242    REAL(rstd) :: temp(iim*jjm,llm)
243    REAL(rstd) :: pks(iim*jjm)
244    REAL(rstd) :: pk(iim*jjm,llm)
245    REAL(rstd) :: ue(3*iim*jjm,llm)
246    REAL(rstd) :: q(iim*jjm,llm,nqtot)
247
248    REAL(rstd) :: p(iim*jjm,llm+1)
249    REAL(rstd) :: uc(iim*jjm,llm,3)
250    REAL(rstd) :: ulon(iim*jjm,llm)
251    REAL(rstd) :: ulat(iim*jjm,llm)
252
253!$OMP BARRIER
254    CALL compute_pression(ps,p,0)
255!$OMP BARRIER
256    CALL compute_exner(ps,p,pks,pk,0) 
257!$OMP BARRIER
258    CALL compute_wind_centered(ue,uc)
259    CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
260!$OMP BARRIER
261    IF (is_omp_level_master) THEN
262      CALL pack_domain(info, phis, physics_inout%phis)
263      CALL pack_domain(info, p, physics_inout%p)
264      CALL pack_domain(info, pk, physics_inout%pk)
265      CALL pack_domain(info, Temp, physics_inout%Temp)
266      CALL pack_domain(info, ulon, physics_inout%ulon)
267      CALL pack_domain(info, ulat, physics_inout%ulat)
268      CALL pack_domain(info, q, physics_inout%q)
269    ENDIF
270!$OMP BARRIER
271  END SUBROUTINE pack_physics
272
273  SUBROUTINE unpack_physics(info, ps,temp, q, dulon, dulat)
274    USE icosa
275    USE physics_interface_mod
276    USE theta2theta_rhodz_mod
277    USE omp_para
278    IMPLICIT NONE
279    TYPE(t_pack_info) :: info
280    REAL(rstd) :: ps(iim*jjm)
281    REAL(rstd) :: temp(iim*jjm,llm)
282    REAL(rstd) :: q(iim*jjm,llm,nqtot)
283    REAL(rstd) :: dulon(iim*jjm,llm)
284    REAL(rstd) :: dulat(iim*jjm,llm)
285
286    REAL(rstd) :: dq(iim*jjm,llm,nqtot)
287    REAL(rstd) :: dTemp(iim*jjm,llm)
288
289!$OMP BARRIER
290    IF (is_omp_level_master) THEN
291      CALL unpack_domain(info, dulon, physics_inout%dulon)
292      CALL unpack_domain(info, dulat, physics_inout%dulat)
293      CALL unpack_domain(info, dq, physics_inout%dq)
294      CALL unpack_domain(info, Temp, physics_inout%Temp)
295      CALL unpack_domain(info, dTemp, physics_inout%dTemp)
296      q = q + physics_inout%dt_phys * dq
297      Temp = Temp + physics_inout%dt_phys * dTemp
298    ENDIF
299!$OMP BARRIER
300
301!    CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
302  END SUBROUTINE unpack_physics
303
304  SUBROUTINE compute_update_velocity(dulon, dulat, ue)
305    USE icosa
306    USE physics_interface_mod
307    USE wind_mod
308    USE omp_para
309    IMPLICIT NONE
310    REAL(rstd) :: dulon(iim*jjm,llm)
311    REAL(rstd) :: dulat(iim*jjm,llm)
312    REAL(rstd) :: ue(3*iim*jjm,llm)
313    REAL(rstd) :: duc(iim*jjm,llm,3)
314    REAL(rstd) :: dt2, due
315    INTEGER :: i,j,ij,l
316    ! Reconstruct wind tendencies at edges and add to normal wind
317    CALL compute_wind_centered_from_lonlat_compound(dulon,dulat,duc)
318    dt2=.5*physics_inout%dt_phys
319    DO l=ll_begin,ll_end
320      DO j=jj_begin,jj_end
321        DO i=ii_begin,ii_end
322          ij=(j-1)*iim+i
323          due = sum( (duc(ij,l,:) + duc(ij+t_right,l,:))*ep_e(ij+u_right,:) )
324          ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due
325
326          due = sum( (duc(ij,l,:) + duc(ij+t_lup,l,:))*ep_e(ij+u_lup,:) )
327          ue(ij+u_lup,l)=ue(ij+u_lup,l) + dt2*due
328
329          due = sum( (duc(ij,l,:) + duc(ij+t_ldown,l,:))*ep_e(ij+u_ldown,:) )
330          ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + dt2*due
331        ENDDO
332      ENDDO
333    ENDDO
334  END SUBROUTINE compute_update_velocity
335
336END MODULE physics_mod
Note: See TracBrowser for help on using the repository browser.