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

Last change on this file since 276 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
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
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    physics_inout%dt_phys = dt*itau_physics
29    physics_type='none'
30    CALL getin("physics",physics_type)
31    SELECT CASE(TRIM(physics_type))
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
37    CASE ('dcmip')
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
41       CALL init_physics_dcmip
42       CALL init_pack_after ! Defines Ai, lon, lat in physics_inout
43       phys_type = phys_DCMIP
44    CASE DEFAULT
45       IF(is_mpi_root) PRINT*, 'init_physics : Bad selector for variable physics <',&
46            TRIM(physics_type), '> options are <none>, <held_suarez>, <dcmip>'
47       STOP
48    END SELECT
49
50    IF(is_mpi_root) PRINT *, 'phys_type = ',phys_type
51  END SUBROUTINE init_physics
52
53  SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
54    USE icosa
55    USE physics_interface_mod
56    USE physics_dcmip_mod, ONLY : write_physics_dcmip => write_physics
57    USE etat0_heldsz_mod
58    IMPLICIT NONE
59    INTEGER, INTENT(IN)   :: it
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(:)
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
71    LOGICAL:: firstcall,lastcall
72    INTEGER :: ind
73    TYPE(t_physics_inout) :: args
74
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
83          CALL physics_column(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
84       END SELECT
85
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
90
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
100  SUBROUTINE physics_column(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
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
134    END SELECT
135
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
147
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
162  END SUBROUTINE physics_column
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
253END MODULE physics_mod
Note: See TracBrowser for help on using the repository browser.