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

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

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