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

Last change on this file since 325 was 325, checked in by dubos, 9 years ago

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