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

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

Merging OpenMP parallisme mode : by subdomain and on vertical level.
This feature is actually experimental but may be retro-compatible with the last method based only on subdomain

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