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

Last change on this file since 275 was 266, checked in by ymipsl, 10 years ago

Synchronize trunk and Saturn branch.
Merge modification from Saturn branch to trunk

YM
File size: 7.8 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
[213]75    IF(MOD(it+1,itau_physics)==0) THEN
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
104    IMPLICIT NONE
105    TYPE(t_field),POINTER :: f_phis(:)
106    TYPE(t_field),POINTER :: f_ps(:)
107    TYPE(t_field),POINTER :: f_theta_rhodz(:)
108    TYPE(t_field),POINTER :: f_ue(:)
109    TYPE(t_field),POINTER :: f_q(:)
110    REAL(rstd),POINTER :: phis(:)
111    REAL(rstd),POINTER :: ps(:)
112    REAL(rstd),POINTER :: theta_rhodz(:,:)
113    REAL(rstd),POINTER :: ue(:,:)
114    REAL(rstd),POINTER :: dulon(:,:)
115    REAL(rstd),POINTER :: dulat(:,:)
116    REAL(rstd),POINTER :: q(:,:,:)
117    INTEGER :: it, ind
118
119    DO ind=1,ndomain
120       IF (.NOT. assigned_domain(ind)) CYCLE
121       CALL swap_dimensions(ind)
122       CALL swap_geometry(ind)
123       phis=f_phis(ind)
124       ps=f_ps(ind)
125       theta_rhodz=f_theta_rhodz(ind)
126       ue=f_ue(ind)
127       q=f_q(ind)
128       CALL pack_physics(pack_info(ind), phis, ps, theta_rhodz, ue, q)
129    END DO
130
131    SELECT CASE(phys_type)
132    CASE (phys_DCMIP)
133       CALL full_physics_dcmip
[196]134    END SELECT
[170]135
[213]136    DO ind=1,ndomain
137       IF (.NOT. assigned_domain(ind)) CYCLE
138       CALL swap_dimensions(ind)
139       CALL swap_geometry(ind)
140       ps=f_ps(ind)
141       theta_rhodz=f_theta_rhodz(ind)
142       q=f_q(ind)
143       dulon=f_dulon(ind)
144       dulat=f_dulat(ind)
145       CALL unpack_physics(pack_info(ind), ps, theta_rhodz, q, dulon, dulat)
146    END DO
[170]147
[213]148    ! Transfer dulon, dulat
149    CALL transfert_request(f_dulon,req_i0)
150    CALL transfert_request(f_dulat,req_i0)
151
152    DO ind=1,ndomain
153       IF (.NOT. assigned_domain(ind)) CYCLE
154       CALL swap_dimensions(ind)
155       CALL swap_geometry(ind)
156       ue=f_ue(ind)
157       dulon=f_dulon(ind)
158       dulat=f_dulat(ind)
159       CALL compute_update_velocity(dulon, dulat, ue)
160    END DO
161
[214]162  END SUBROUTINE physics_column
[213]163
164  SUBROUTINE pack_physics(info, phis, ps, theta_rhodz, ue, q)
165    USE icosa
166    USE wind_mod
167    USE pression_mod
168    USE theta2theta_rhodz_mod
169    USE physics_interface_mod
170    IMPLICIT NONE
171    TYPE(t_pack_info) :: info
172    REAL(rstd) :: phis(iim*jjm)
173    REAL(rstd) :: ps(iim*jjm)
174    REAL(rstd) :: theta_rhodz(iim*jjm,llm)
175    REAL(rstd) :: ue(3*iim*jjm,llm)
176    REAL(rstd) :: q(iim*jjm,llm,nqtot)
177
178    REAL(rstd) :: p(iim*jjm,llm+1)
179    REAL(rstd) :: Temp(iim*jjm,llm)
180    REAL(rstd) :: uc(iim*jjm,3,llm)
181    REAL(rstd) :: ulon(iim*jjm,llm)
182    REAL(rstd) :: ulat(iim*jjm,llm)
183
184    CALL compute_pression(ps,p,0)
185    CALL compute_theta_rhodz2temperature(ps,theta_rhodz,Temp,0)
186    CALL compute_wind_centered(ue,uc)
187    CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
188
189    CALL pack_domain(info, phis, physics_inout%phis)
190    CALL pack_domain(info, p, physics_inout%p)
191    CALL pack_domain(info, Temp, physics_inout%Temp)
192    CALL pack_domain(info, ulon, physics_inout%ulon)
193    CALL pack_domain(info, ulat, physics_inout%ulat)
194    CALL pack_domain(info, q, physics_inout%q)
195  END SUBROUTINE pack_physics
196
197  SUBROUTINE unpack_physics(info, ps,theta_rhodz, q, dulon, dulat)
198    USE icosa
199    USE physics_interface_mod
200    USE theta2theta_rhodz_mod
201    IMPLICIT NONE
202    TYPE(t_pack_info) :: info
203    REAL(rstd) :: ps(iim*jjm)
204    REAL(rstd) :: theta_rhodz(iim*jjm,llm)
205    REAL(rstd) :: Temp(iim*jjm,llm)
206    REAL(rstd) :: q(iim*jjm,llm,nqtot)
207    REAL(rstd) :: dulon(iim*jjm,llm)
208    REAL(rstd) :: dulat(iim*jjm,llm)
209
210    REAL(rstd) :: dq(iim*jjm,llm,nqtot)
211    REAL(rstd) :: dTemp(iim*jjm,llm)
212    CALL unpack_domain(info, dulon, physics_inout%dulon)
213    CALL unpack_domain(info, dulat, physics_inout%dulat)
214    CALL unpack_domain(info, dq, physics_inout%dq)
215    CALL unpack_domain(info, Temp, physics_inout%Temp)
216    CALL unpack_domain(info, dTemp, physics_inout%dTemp)
217    q = q + physics_inout%dt_phys * dq
218    Temp = Temp + physics_inout%dt_phys * dTemp
219    CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
220  END SUBROUTINE unpack_physics
221
222  SUBROUTINE compute_update_velocity(dulon, dulat, ue)
223    USE icosa
224    USE physics_interface_mod
225    USE wind_mod
226    IMPLICIT NONE
227    REAL(rstd) :: dulon(iim*jjm,llm)
228    REAL(rstd) :: dulat(iim*jjm,llm)
229    REAL(rstd) :: ue(3*iim*jjm,llm)
230    REAL(rstd) :: duc(iim*jjm,3,llm)
231    REAL(rstd) :: dt2, due
232    INTEGER :: i,j,ij,l
233    ! Reconstruct wind tendencies at edges and add to normal wind
234    CALL compute_wind_centered_from_lonlat_compound(dulon,dulat,duc)
235    dt2=.5*physics_inout%dt_phys
236    DO l=1,llm
237      DO j=jj_begin,jj_end
238        DO i=ii_begin,ii_end
239          ij=(j-1)*iim+i
240          due = sum( (duc(ij,:,l) + duc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
241          ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due
242
243          due = sum( (duc(ij,:,l) + duc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )
244          ue(ij+u_lup,l)=ue(ij+u_lup,l) + dt2*due
245
246          due = sum( (duc(ij,:,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )
247          ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + dt2*due
248        ENDDO
249      ENDDO
250    ENDDO
251  END SUBROUTINE compute_update_velocity
252
[81]253END MODULE physics_mod
Note: See TracBrowser for help on using the repository browser.