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

Last change on this file since 327 was 327, checked in by ymipsl, 9 years ago

Merge recent developments from saturn branch onto trunk.

  • lmdz generic physics interface
  • performance improvment on mix mpi/openmp
  • asynchrone and overlaping communication
  • best domain distribution between process and threads
  • ....

This version is compatible with the actual saturn version and the both branches are considered merged on dynamico component.

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