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

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

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