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

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

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