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

Last change on this file since 213 was 213, checked in by dubos, 10 years ago

New dyn/phys interface - halo points not passed to physics any more (cleanup follows)
File size: 12.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
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
13  CHARACTER(LEN=255) :: physics_type="automatic"
14!$OMP THREADPRIVATE(physics_type)
15
16  PUBLIC :: physics, init_physics
17
18CONTAINS
19
20  SUBROUTINE init_physics
21    USE mpipara
22    USE etat0_mod
23    USE icosa
24    USE physics_interface_mod
25    USE physics_dcmip_mod, ONLY : &
26         init_physics_dcmip=>init_physics, init_physics_dcmip_new=>init_physics_new
27    IMPLICIT NONE
28
29    CALL allocate_field(f_dulon,field_t,type_real,llm, name='dulon')
30    CALL allocate_field(f_dulat,field_t,type_real,llm, name='dulat')
31    CALL init_pack_before ! Compute physics_inout%ngrid and offsets used by pack/unpack
32   
33    physics_type='automatic'
34    CALL getin("physics",physics_type)
35
36    SELECT CASE(TRIM(physics_type))
37    CASE ('automatic')
38       etat0_type='jablonowsky06'
39       CALL getin("etat0",etat0_type)
40       SELECT CASE(TRIM(etat0_type))
41       CASE('held_suarez')
42          phys_type = phys_HS94
43       CASE DEFAULT
44          IF(is_mpi_root) PRINT*,"NO PHYSICAL PACKAGE USED"
45          phys_type = phys_none
46       END SELECT
47
48    CASE ('dcmip')
49       CALL init_physics_dcmip_new
50!       CALL init_physics_dcmip
51       phys_type = phys_DCMIP
52    CASE DEFAULT
53       IF(is_mpi_root) PRINT*, 'init_physics : Bad selector for variable physics <',&
54            TRIM(physics_type), '> options are <automatic>, <dcmip>, <dry>'
55       STOP
56    END SELECT
57
58    IF(is_mpi_root) THEN
59       PRINT *, 'phys_type = ',phys_type
60       PRINT *, 'nb_extra_physics_2D = ', nb_extra_physics_2D
61       PRINT *, 'nb_extra_physics_3D = ', nb_extra_physics_3D
62    END IF
63
64    IF(.FALSE.) THEN ! draft interface
65       IF(nb_extra_physics_2D>0) CALL allocate_field(f_extra_physics_2D,field_t,type_real,nb_extra_physics_2D)
66       IF(nb_extra_physics_3D>0) CALL allocate_field(f_extra_physics_3D,field_t,type_real,llm,nb_extra_physics_3D)
67    ELSE
68       CALL init_pack_after ! Defines Ai, lon, lat in physics_inout
69    END IF
70  END SUBROUTINE init_physics
71
72  SUBROUTINE physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
73    USE icosa
74    USE physics_interface_mod
75    USE physics_dcmip_mod, ONLY : write_physics_dcmip => write_physics
76    USE etat0_heldsz_mod
77    IMPLICIT NONE
78    INTEGER, INTENT(IN)   :: it
79    REAL(rstd),INTENT(IN)::jD_cur,jH_cur
80    TYPE(t_field),POINTER :: f_phis(:)
81    TYPE(t_field),POINTER :: f_ps(:)
82    TYPE(t_field),POINTER :: f_theta_rhodz(:)
83    TYPE(t_field),POINTER :: f_ue(:)
84    TYPE(t_field),POINTER :: f_q(:)
85    REAL(rstd),POINTER :: phis(:)
86    REAL(rstd),POINTER :: ps(:)
87    REAL(rstd),POINTER :: theta_rhodz(:,:)
88    REAL(rstd),POINTER :: ue(:,:)
89    REAL(rstd),POINTER :: q(:,:,:)
90
91    LOGICAL:: firstcall,lastcall
92    INTEGER :: ind
93    TYPE(t_physics_inout) :: args
94
95    IF(MOD(it+1,itau_physics)==0) THEN
96   
97       SELECT CASE(phys_type)
98       CASE (phys_none)
99          ! No physics, do nothing
100       CASE(phys_HS94)
101          CALL held_suarez(f_ps,f_theta_rhodz,f_ue) 
102       CASE DEFAULT
103          IF(.FALSE.) THEN ! draft interface
104             args%dt_phys = dt*itau_physics
105             DO ind=1,ndomain
106                IF (.NOT. assigned_domain(ind)) CYCLE
107                CALL swap_dimensions(ind)
108                CALL swap_geometry(ind)
109               
110                phis=f_phis(ind)
111                ps=f_ps(ind)
112                theta_rhodz=f_theta_rhodz(ind)
113                ue=f_ue(ind)
114                q=f_q(ind)
115               
116                IF(nb_extra_physics_2D>0) args%extra_2D=f_extra_physics_2D(ind)
117                IF(nb_extra_physics_3D>0) args%extra_3D=f_extra_physics_3D(ind)
118                CALL physics_column(args, phis, ps, theta_rhodz, ue, q)
119             ENDDO
120             
121             IF (mod(it,itau_out)==0 ) THEN
122                IF(nb_extra_physics_2D>0) CALL writefield("extra_physics_2D",f_extra_physics_2D)
123                IF(nb_extra_physics_3D>0) CALL writefield("extra_physics_3D",f_extra_physics_3D)
124             ENDIF
125          ELSE ! new interface
126             physics_inout%dt_phys = dt*itau_physics
127             CALL physics_column_new(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
128          END IF
129       END SELECT
130
131       CALL transfert_request(f_theta_rhodz,req_i0)
132       CALL transfert_request(f_ue,req_e0_vect)
133       CALL transfert_request(f_q,req_i0)
134    END IF
135
136    IF (mod(it,itau_out)==0 ) THEN
137       SELECT CASE(phys_type)
138       CASE (phys_DCMIP)
139          CALL write_physics_dcmip
140       END SELECT
141    END IF
142   
143  END SUBROUTINE physics
144
145!--------------------------------- New interface --------------------------------------
146
147  SUBROUTINE physics_column_new(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
148    USE icosa
149    USE physics_interface_mod
150    USE physics_dcmip_mod, ONLY : full_physics_dcmip => full_physics
151    IMPLICIT NONE
152    TYPE(t_field),POINTER :: f_phis(:)
153    TYPE(t_field),POINTER :: f_ps(:)
154    TYPE(t_field),POINTER :: f_theta_rhodz(:)
155    TYPE(t_field),POINTER :: f_ue(:)
156    TYPE(t_field),POINTER :: f_q(:)
157    REAL(rstd),POINTER :: phis(:)
158    REAL(rstd),POINTER :: ps(:)
159    REAL(rstd),POINTER :: theta_rhodz(:,:)
160    REAL(rstd),POINTER :: ue(:,:)
161    REAL(rstd),POINTER :: dulon(:,:)
162    REAL(rstd),POINTER :: dulat(:,:)
163    REAL(rstd),POINTER :: q(:,:,:)
164    INTEGER :: it, ind
165
166    DO ind=1,ndomain
167       IF (.NOT. assigned_domain(ind)) CYCLE
168       CALL swap_dimensions(ind)
169       CALL swap_geometry(ind)
170       phis=f_phis(ind)
171       ps=f_ps(ind)
172       theta_rhodz=f_theta_rhodz(ind)
173       ue=f_ue(ind)
174       q=f_q(ind)
175       CALL pack_physics(pack_info(ind), phis, ps, theta_rhodz, ue, q)
176    END DO
177
178    SELECT CASE(phys_type)
179    CASE (phys_DCMIP)
180       CALL full_physics_dcmip
181    END SELECT
182
183    DO ind=1,ndomain
184       IF (.NOT. assigned_domain(ind)) CYCLE
185       CALL swap_dimensions(ind)
186       CALL swap_geometry(ind)
187       ps=f_ps(ind)
188       theta_rhodz=f_theta_rhodz(ind)
189       q=f_q(ind)
190       dulon=f_dulon(ind)
191       dulat=f_dulat(ind)
192       CALL unpack_physics(pack_info(ind), ps, theta_rhodz, q, dulon, dulat)
193    END DO
194
195    ! Transfer dulon, dulat
196    CALL transfert_request(f_dulon,req_i0)
197    CALL transfert_request(f_dulat,req_i0)
198
199    DO ind=1,ndomain
200       IF (.NOT. assigned_domain(ind)) CYCLE
201       CALL swap_dimensions(ind)
202       CALL swap_geometry(ind)
203       ue=f_ue(ind)
204       dulon=f_dulon(ind)
205       dulat=f_dulat(ind)
206       CALL compute_update_velocity(dulon, dulat, ue)
207    END DO
208
209  END SUBROUTINE physics_column_new
210
211  SUBROUTINE pack_physics(info, phis, ps, theta_rhodz, ue, q)
212    USE icosa
213    USE wind_mod
214    USE pression_mod
215    USE theta2theta_rhodz_mod
216    USE physics_interface_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) :: theta_rhodz(iim*jjm,llm)
222    REAL(rstd) :: ue(3*iim*jjm,llm)
223    REAL(rstd) :: q(iim*jjm,llm,nqtot)
224
225    REAL(rstd) :: p(iim*jjm,llm+1)
226    REAL(rstd) :: Temp(iim*jjm,llm)
227    REAL(rstd) :: uc(iim*jjm,3,llm)
228    REAL(rstd) :: ulon(iim*jjm,llm)
229    REAL(rstd) :: ulat(iim*jjm,llm)
230
231    CALL compute_pression(ps,p,0)
232    CALL compute_theta_rhodz2temperature(ps,theta_rhodz,Temp,0)
233    CALL compute_wind_centered(ue,uc)
234    CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
235
236    CALL pack_domain(info, phis, physics_inout%phis)
237    CALL pack_domain(info, p, physics_inout%p)
238    CALL pack_domain(info, Temp, physics_inout%Temp)
239    CALL pack_domain(info, ulon, physics_inout%ulon)
240    CALL pack_domain(info, ulat, physics_inout%ulat)
241    CALL pack_domain(info, q, physics_inout%q)
242  END SUBROUTINE pack_physics
243
244  SUBROUTINE unpack_physics(info, ps,theta_rhodz, q, dulon, dulat)
245    USE icosa
246    USE physics_interface_mod
247    USE theta2theta_rhodz_mod
248    IMPLICIT NONE
249    TYPE(t_pack_info) :: info
250    REAL(rstd) :: ps(iim*jjm)
251    REAL(rstd) :: theta_rhodz(iim*jjm,llm)
252    REAL(rstd) :: Temp(iim*jjm,llm)
253    REAL(rstd) :: q(iim*jjm,llm,nqtot)
254    REAL(rstd) :: dulon(iim*jjm,llm)
255    REAL(rstd) :: dulat(iim*jjm,llm)
256
257    REAL(rstd) :: dq(iim*jjm,llm,nqtot)
258    REAL(rstd) :: dTemp(iim*jjm,llm)
259    CALL unpack_domain(info, dulon, physics_inout%dulon)
260    CALL unpack_domain(info, dulat, physics_inout%dulat)
261    CALL unpack_domain(info, dq, physics_inout%dq)
262    CALL unpack_domain(info, Temp, physics_inout%Temp)
263    CALL unpack_domain(info, dTemp, physics_inout%dTemp)
264    q = q + physics_inout%dt_phys * dq
265    Temp = Temp + physics_inout%dt_phys * dTemp
266    CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
267  END SUBROUTINE unpack_physics
268
269  SUBROUTINE compute_update_velocity(dulon, dulat, ue)
270    USE icosa
271    USE physics_interface_mod
272    USE wind_mod
273    IMPLICIT NONE
274    REAL(rstd) :: dulon(iim*jjm,llm)
275    REAL(rstd) :: dulat(iim*jjm,llm)
276    REAL(rstd) :: ue(3*iim*jjm,llm)
277    REAL(rstd) :: duc(iim*jjm,3,llm)
278    REAL(rstd) :: dt2, due
279    INTEGER :: i,j,ij,l
280    ! Reconstruct wind tendencies at edges and add to normal wind
281    CALL compute_wind_centered_from_lonlat_compound(dulon,dulat,duc)
282    dt2=.5*physics_inout%dt_phys
283    DO l=1,llm
284      DO j=jj_begin,jj_end
285        DO i=ii_begin,ii_end
286          ij=(j-1)*iim+i
287          due = sum( (duc(ij,:,l) + duc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
288          ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due
289
290          due = sum( (duc(ij,:,l) + duc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )
291          ue(ij+u_lup,l)=ue(ij+u_lup,l) + dt2*due
292
293          due = sum( (duc(ij,:,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )
294          ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + dt2*due
295        ENDDO
296      ENDDO
297    ENDDO
298  END SUBROUTINE compute_update_velocity
299
300!--------------------------------- Draft interface --------------------------------------
301
302  SUBROUTINE physics_column(args, phis, ps, theta_rhodz, ue, q)
303    USE icosa
304    USE wind_mod
305    USE pression_mod
306    USE theta2theta_rhodz_mod
307    USE physics_interface_mod
308    USE physics_dcmip_mod
309    IMPLICIT NONE
310    TYPE(t_physics_inout) :: args   
311    REAL(rstd) :: phis(iim*jjm)
312    REAL(rstd) :: ps(iim*jjm)
313    REAL(rstd) :: theta_rhodz(iim*jjm,llm)
314    REAL(rstd) :: ue(3*iim*jjm,llm)
315    REAL(rstd), TARGET :: q(iim*jjm,llm,nqtot)
316    ! local arrays
317    REAL(rstd), TARGET :: lat(iim*jjm)
318    REAL(rstd), TARGET :: lon(iim*jjm)
319    REAL(rstd), TARGET :: p(iim*jjm,llm+1)
320    REAL(rstd), TARGET :: Temp(iim*jjm,llm)
321    REAL(rstd), TARGET :: ulon(iim*jjm,llm)
322    REAL(rstd), TARGET :: ulat(iim*jjm,llm)
323    REAL(rstd), TARGET :: dTemp(iim*jjm,llm)
324    REAL(rstd), TARGET :: dulon(iim*jjm,llm)
325    REAL(rstd), TARGET :: dulat(iim*jjm,llm)
326    REAL(rstd), TARGET :: dq(iim*jjm,llm,nqtot)
327    REAL(rstd) :: uc(iim*jjm,3,llm)  ! 3D velocity at cell centers
328
329    INTEGER :: i,j,ij,l
330    REAL(rstd) :: due, dt2
331
332    DO j=jj_begin,jj_end
333      DO i=ii_begin,ii_end
334        ij=(j-1)*iim+i
335        CALL xyz2lonlat(xyz_i(ij,:),lon(ij),lat(ij)) 
336      ENDDO
337    ENDDO
338
339    ! Reconstruct wind vector at hexagons
340    CALL compute_pression(ps,p,0)
341    CALL compute_theta_rhodz2temperature(ps,theta_rhodz,Temp,0)
342    CALL compute_wind_centered(ue,uc)
343    CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
344    args%ngrid = iim*jjm
345    args%lon => lon
346    args%lat => lat
347    args%p => p
348    args%Temp => Temp
349    args%ulon => ulon
350    args%ulat => ulat
351    args%q => q
352    args%dTemp => dTemp
353    args%dulon => dulon
354    args%dulat => dulat
355    args%dq => dq
356
357    SELECT CASE(phys_type)
358    CASE (phys_DCMIP)
359       CALL compute_phys_wrap(args)
360    END SELECT
361
362    q = q + args%dt_phys * dq
363    Temp = Temp + args%dt_phys * dTemp
364    CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
365   
366    ! Reconstruct wind tendencies at edges and add
367    CALL compute_wind_centered_from_lonlat_compound(dulon,dulat,uc)
368    dt2=.5*args%dt_phys
369    DO l=1,llm
370      DO j=jj_begin,jj_end
371        DO i=ii_begin,ii_end
372          ij=(j-1)*iim+i
373          due = sum( (uc(ij,:,l) + uc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
374          ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due
375
376          due = sum( (uc(ij,:,l) + uc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )
377          ue(ij+u_lup,l)=ue(ij+u_lup,l) + dt2*due
378
379          due = sum( (uc(ij,:,l) + uc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )
380          ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + dt2*due
381        ENDDO
382      ENDDO
383    ENDDO
384
385  END SUBROUTINE physics_column
386
387END MODULE physics_mod
Note: See TracBrowser for help on using the repository browser.