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

Last change on this file since 899 was 899, checked in by adurocher, 5 years ago

trunk : Fixed GCC warnings

Fixed iso c bindings
fixed warnings with -Wall -Wno-aliasing -Wno-unused -Wno-unused-dummy-argument -Wno-maybe-uninitialized -Wno-tabs warnings
Removed all unused variables (-Wunused-variable)
vector%dot_product is now dot_product_3d to avoid compilation warning "dot_product shadows intrinsic" with GCC
File size: 13.1 KB
Line 
1MODULE physics_mod
2  USE icosa
3  USE field_mod
4  USE physics_interface_mod
5  USE omp_para
6  IMPLICIT NONE
7  PRIVATE
8
9  INTEGER, PARAMETER :: phys_none=0, phys_HS94=1, phys_DCMIP=2, phys_lmdz_generic=3, phys_LB2012=4, phys_external=5,&
10                        phys_DCMIP2016=6
11
12  INTEGER :: phys_type
13  TYPE(t_field),POINTER,SAVE :: f_dulon(:), f_dulat(:)
14  TYPE(t_field),POINTER,SAVE :: f_ulon(:), f_ulat(:)
15  TYPE(t_field),POINTER,SAVE :: f_p(:), f_pk(:)
16  TYPE(t_field),POINTER,SAVE :: f_temp(:)
17  TYPE(t_field),POINTER,SAVE :: f_du_phys(:)
18
19  CHARACTER(LEN=255),SAVE :: physics_type
20!$OMP THREADPRIVATE(physics_type)
21  TYPE(t_message),SAVE :: req_theta0, req_ue0, req_q0
22
23  PUBLIC :: physics, init_physics, zero_du_phys
24
25CONTAINS
26
27  SUBROUTINE init_physics
28    USE mpipara
29    USE etat0_mod
30    USE physics_dcmip_mod, ONLY : init_physics_dcmip=>init_physics
31    USE physics_dcmip2016_mod, ONLY : init_physics_dcmip2016=>init_physics
32    USE etat0_venus_mod, ONLY : init_phys_venus=>init_physics
33    USE physics_lmdz_generic_mod, ONLY : init_physics_lmdz_generic=>init_physics
34    USE physics_external_mod, ONLY : init_physics_external=>init_physics
35
36    physics_inout%dt_phys = dt*itau_physics
37    physics_type='none'
38    CALL getin("physics",physics_type)
39    SELECT CASE(TRIM(physics_type))
40    CASE ('none')
41
42!$OMP PARALLEL
43       IF(is_mpi_root) PRINT*,"NO PHYSICAL PACKAGE USED"
44       phys_type = phys_none
45!$OMP END PARALLEL
46
47    CASE ('held_suarez')
48
49!$OMP PARALLEL
50       phys_type = phys_HS94
51!$OMP END PARALLEL
52
53    CASE ('Lebonnois2012')
54
55!$OMP PARALLEL
56       phys_type = phys_LB2012
57       CALL init_phys_venus
58!$OMP END PARALLEL
59
60    CASE ('phys_lmdz_generic')
61
62!$OMP PARALLEL
63       CALL init_physics_lmdz_generic
64       phys_type=phys_lmdz_generic
65!$OMP END PARALLEL
66
67    CASE ('phys_external')
68
69       CALL init_physics_external
70!$OMP PARALLEL
71       phys_type=phys_external
72!$OMP END PARALLEL
73
74    CASE ('dcmip')
75
76!$OMP PARALLEL
77       CALL allocate_field(f_dulon,field_t,type_real,llm, name='dulon')
78       CALL allocate_field(f_dulat,field_t,type_real,llm, name='dulat')
79       CALL allocate_field(f_temp,field_t,type_real,llm, name='temp')
80       CALL allocate_field(f_ulon,field_t,type_real,llm, name='ulon')
81       CALL allocate_field(f_ulat,field_t,type_real,llm, name='ulat')
82       CALL allocate_field(f_p,field_t,type_real,llm+1, name='p')
83       CALL allocate_field(f_pk,field_t,type_real,llm, name='pk')
84       CALL init_pack_before ! Compute physics_inout%ngrid and offsets used by pack/unpack
85       CALL init_physics_dcmip
86       CALL init_pack_after ! Defines Ai, lon, lat in physics_inout
87       phys_type = phys_DCMIP
88!$OMP END PARALLEL
89
90    CASE ('dcmip2016')
91
92!$OMP PARALLEL
93       CALL allocate_field(f_dulon,field_t,type_real,llm, name='dulon')
94       CALL allocate_field(f_dulat,field_t,type_real,llm, name='dulat')
95       CALL allocate_field(f_temp,field_t,type_real,llm, name='temp')
96       CALL allocate_field(f_ulon,field_t,type_real,llm, name='ulon')
97       CALL allocate_field(f_ulat,field_t,type_real,llm, name='ulat')
98       CALL allocate_field(f_p,field_t,type_real,llm+1, name='p')
99       CALL allocate_field(f_pk,field_t,type_real,llm, name='pk')
100       CALL init_pack_before ! Compute physics_inout%ngrid and offsets used by pack/unpack
101       CALL init_physics_dcmip2016
102       CALL init_pack_after ! Defines Ai, lon, lat in physics_inout
103       phys_type = phys_DCMIP2016
104!$OMP END PARALLEL
105
106    CASE DEFAULT
107       IF(is_mpi_root) PRINT*, 'init_physics : Bad selector for variable physics <',&
108            TRIM(physics_type), '> options are <none>, <held_suarez>, <Lebonnois2012>, <dcmip>', &
109                                '<phys_lmdz_generic>, <phys_external>'
110       STOP
111    END SELECT
112
113!$OMP PARALLEL
114    CALL allocate_field(f_du_phys,field_u,type_real,llm, name='du_phys')
115
116    IF(is_mpi_root) PRINT *, 'phys_type = ',phys_type
117   
118   
119!$OMP END PARALLEL
120  END SUBROUTINE init_physics
121
122  SUBROUTINE zero_du_phys()
123    REAL(rstd), DIMENSION(:,:), POINTER :: du
124    INTEGER :: ind
125    DO ind=1,ndomain
126       IF (.NOT. assigned_domain(ind)) CYCLE
127       CALL swap_dimensions(ind)
128       CALL swap_geometry(ind)
129       du=f_du_phys(ind)
130       du(:,ll_begin:ll_end) = 0.
131    END DO
132  END SUBROUTINE zero_du_phys
133
134  SUBROUTINE add_du_phys(coef, f_u)
135    REAL(rstd), INTENT(IN) :: coef  ! -1 before physics, +1 after physics
136    TYPE(t_field),POINTER :: f_u(:) ! velocity field before/after call to physics
137    REAL(rstd), DIMENSION(:,:), POINTER :: u, du
138    INTEGER :: ind
139    DO ind=1,ndomain
140       IF (.NOT. assigned_domain(ind)) CYCLE
141       CALL swap_dimensions(ind)
142       CALL swap_geometry(ind)
143       du=f_du_phys(ind)
144       u=f_u(ind)
145       du(:,ll_begin:ll_end) = du(:,ll_begin:ll_end) + coef*u(:,ll_begin:ll_end)
146    END DO
147  END SUBROUTINE add_du_phys
148
149  SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
150    USE physics_lmdz_generic_mod, ONLY : physics_lmdz_generic => physics
151    USE physics_external_mod, ONLY : physics_external => physics
152    USE physics_dcmip_mod, ONLY : write_physics_dcmip => write_physics
153    USE physics_dcmip2016_mod, ONLY : write_physics_dcmip2016 => write_physics
154    USE etat0_heldsz_mod
155    USE etat0_venus_mod, ONLY : phys_venus => physics
156    INTEGER, INTENT(IN)   :: it
157    TYPE(t_field),POINTER :: f_phis(:)
158    TYPE(t_field),POINTER :: f_ps(:)
159    TYPE(t_field),POINTER :: f_theta_rhodz(:)
160    TYPE(t_field),POINTER :: f_ue(:)
161    TYPE(t_field),POINTER :: f_wflux(:)
162    TYPE(t_field),POINTER :: f_q(:)
163
164    LOGICAL,SAVE :: first=.TRUE.
165    !$OMP THREADPRIVATE(first)
166    IF (first) THEN
167      CALL init_message(f_theta_rhodz, req_i0, req_theta0)
168      CALL init_message(f_ue, req_e0_vect, req_ue0)
169      CALL init_message(f_q, req_i0, req_q0)
170      first=.FALSE.
171    ENDIF
172   
173    IF (phys_external) THEN
174   
175      CALL physics_external(it ,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
176
177    ELSE
178   
179      IF(MOD(it,itau_physics)==0 .AND. phys_type/=phys_none) THEN
180
181         ! as a result of the the two calls to add_du_phys,
182         ! du_phys increases by u(after physics) - u (before physics)
183         CALL add_du_phys(-1., f_ue)
184
185         SELECT CASE(phys_type)
186         CASE(phys_HS94)
187            CALL held_suarez(f_ps,f_theta_rhodz,f_ue) 
188         CASE (phys_lmdz_generic)
189           CALL physics_lmdz_generic(it ,f_phis, f_ps, f_theta_rhodz, f_ue, f_wflux, f_q)
190         CASE(phys_LB2012)
191            CALL phys_venus(f_ps,f_theta_rhodz,f_ue) 
192         CASE DEFAULT
193            CALL physics_column(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
194         END SELECT
195
196         CALL send_message(f_theta_rhodz, req_theta0)
197         CALL send_message(f_ue, req_ue0)
198         CALL send_message(f_q, req_q0)
199         CALL wait_message(req_theta0)
200         CALL wait_message(req_ue0)
201         CALL wait_message(req_q0)
202       
203         CALL add_du_phys(1., f_ue)
204      END IF
205
206      IF (mod(it,itau_out)==0 ) THEN
207         CALL write_physics_tendencies
208         CALL zero_du_phys
209         SELECT CASE(phys_type)
210         CASE (phys_DCMIP)
211            CALL write_physics_dcmip
212         CASE (phys_DCMIP2016)
213            CALL write_physics_dcmip2016
214         END SELECT
215      END IF
216    ENDIF
217   
218  END SUBROUTINE physics
219
220  SUBROUTINE write_physics_tendencies
221    USE observable_mod, ONLY : f_buf_ulon, f_buf_ulat
222    USE wind_mod
223    USE output_field_mod
224    CALL transfert_request(f_du_phys,req_e1_vect)
225    CALL un2ulonlat(f_du_phys, f_buf_ulon, f_buf_ulat, (1./(dt*itau_out)))
226    CALL output_field("dulon_phys",f_buf_ulon)
227    CALL output_field("dulat_phys",f_buf_ulat)
228  END SUBROUTINE write_physics_tendencies
229   
230  SUBROUTINE physics_column(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
231    USE physics_dcmip_mod, ONLY : full_physics_dcmip => full_physics
232    USE physics_dcmip2016_mod, ONLY : full_physics_dcmip2016 => full_physics
233    USE theta2theta_rhodz_mod
234    USE mpipara
235    USE checksum_mod
236    TYPE(t_field),POINTER :: f_phis(:)
237    TYPE(t_field),POINTER :: f_ps(:)
238    TYPE(t_field),POINTER :: f_theta_rhodz(:)
239    TYPE(t_field),POINTER :: f_ue(:)
240    TYPE(t_field),POINTER :: f_q(:)
241    REAL(rstd),POINTER :: phis(:)
242    REAL(rstd),POINTER :: ps(:)
243    REAL(rstd),POINTER :: temp(:,:)
244    REAL(rstd),POINTER :: ue(:,:)
245    REAL(rstd),POINTER :: dulon(:,:)
246    REAL(rstd),POINTER :: dulat(:,:)
247    REAL(rstd),POINTER :: q(:,:,:)
248    REAL(rstd),POINTER :: p(:,:)
249    REAL(rstd),POINTER :: pk(:,:)
250    REAL(rstd),POINTER :: ulon(:,:)
251    REAL(rstd),POINTER :: ulat(:,:)
252    INTEGER :: it, ind
253
254    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_temp)
255   
256    DO ind=1,ndomain
257       IF (.NOT. assigned_domain(ind)) CYCLE
258       CALL swap_dimensions(ind)
259       CALL swap_geometry(ind)
260       phis=f_phis(ind)
261       ps=f_ps(ind)
262       temp=f_temp(ind)
263       ue=f_ue(ind)
264       q=f_q(ind)
265       p=f_p(ind)
266       pk=f_pk(ind)
267       ulon=f_ulon(ind)
268       ulat=f_ulat(ind)
269       CALL pack_physics(pack_info(ind), phis, ps, temp, ue, q, p, pk, ulon, ulat)
270    END DO
271
272    SELECT CASE(phys_type)
273    CASE (phys_DCMIP)
274       IF (is_omp_level_master) CALL full_physics_dcmip
275    CASE (phys_DCMIP2016)
276       IF (is_omp_level_master) CALL full_physics_dcmip2016
277    CASE DEFAULT
278       IF(is_master) PRINT *,'Internal error : illegal value of phys_type', phys_type
279       STOP
280    END SELECT
281
282    DO ind=1,ndomain
283       IF (.NOT. assigned_domain(ind)) CYCLE
284       CALL swap_dimensions(ind)
285       CALL swap_geometry(ind)
286       ps=f_ps(ind)
287       temp=f_temp(ind)
288       q=f_q(ind)
289       dulon=f_dulon(ind)
290       dulat=f_dulat(ind)
291       CALL unpack_physics(pack_info(ind), ps, temp, q, dulon, dulat)
292    END DO
293   
294    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
295
296    ! Transfer dulon, dulat
297    CALL transfert_request(f_dulon,req_i0)
298    CALL transfert_request(f_dulat,req_i0)
299
300    DO ind=1,ndomain
301       IF (.NOT. assigned_domain(ind)) CYCLE
302       CALL swap_dimensions(ind)
303       CALL swap_geometry(ind)
304       ue=f_ue(ind)
305       dulon=f_dulon(ind)
306       dulat=f_dulat(ind)
307       CALL compute_update_velocity(dulon, dulat, ue)
308    END DO
309
310  END SUBROUTINE physics_column
311
312  SUBROUTINE pack_physics(info, phis, ps, temp, ue, q, p, pk, ulon, ulat )
313    USE wind_mod
314    USE pression_mod
315    USE theta2theta_rhodz_mod
316    USE exner_mod
317    TYPE(t_pack_info) :: info
318    REAL(rstd) :: phis(iim*jjm)
319    REAL(rstd) :: ps(iim*jjm)
320    REAL(rstd) :: temp(iim*jjm,llm)
321    REAL(rstd) :: pks(iim*jjm)
322    REAL(rstd) :: pk(iim*jjm,llm)
323    REAL(rstd) :: ue(3*iim*jjm,llm)
324    REAL(rstd) :: q(iim*jjm,llm,nqtot)
325
326    REAL(rstd) :: p(iim*jjm,llm+1)
327    REAL(rstd) :: uc(iim*jjm,llm,3)
328    REAL(rstd) :: ulon(iim*jjm,llm)
329    REAL(rstd) :: ulat(iim*jjm,llm)
330
331!$OMP BARRIER
332    CALL compute_pression(ps,p,0)
333!$OMP BARRIER
334    CALL compute_exner(ps,p,pks,pk,0) 
335!$OMP BARRIER
336    CALL compute_wind_centered(ue,uc)
337    CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
338!$OMP BARRIER
339    IF (is_omp_level_master) THEN
340      CALL pack_domain(info, phis, physics_inout%phis)
341      CALL pack_domain(info, p, physics_inout%p)
342      CALL pack_domain(info, pk, physics_inout%pk)
343      CALL pack_domain(info, Temp, physics_inout%Temp)
344      CALL pack_domain(info, ulon, physics_inout%ulon)
345      CALL pack_domain(info, ulat, physics_inout%ulat)
346      CALL pack_domain(info, q, physics_inout%q)
347    ENDIF
348!$OMP BARRIER
349  END SUBROUTINE pack_physics
350
351  SUBROUTINE unpack_physics(info, ps,temp, q, dulon, dulat)
352    USE theta2theta_rhodz_mod
353    TYPE(t_pack_info) :: info
354    REAL(rstd) :: ps(iim*jjm)
355    REAL(rstd) :: temp(iim*jjm,llm)
356    REAL(rstd) :: q(iim*jjm,llm,nqtot)
357    REAL(rstd) :: dulon(iim*jjm,llm)
358    REAL(rstd) :: dulat(iim*jjm,llm)
359
360    REAL(rstd) :: dq(iim*jjm,llm,nqtot)
361    REAL(rstd) :: dTemp(iim*jjm,llm)
362
363!$OMP BARRIER
364    IF (is_omp_level_master) THEN
365      CALL unpack_domain(info, dulon, physics_inout%dulon)
366      CALL unpack_domain(info, dulat, physics_inout%dulat)
367      CALL unpack_domain(info, dq, physics_inout%dq)
368      CALL unpack_domain(info, Temp, physics_inout%Temp)
369      CALL unpack_domain(info, dTemp, physics_inout%dTemp)
370      q = q + physics_inout%dt_phys * dq
371      Temp = Temp + physics_inout%dt_phys * dTemp
372    ENDIF
373!$OMP BARRIER
374
375!    CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
376  END SUBROUTINE unpack_physics
377
378  SUBROUTINE compute_update_velocity(dulon, dulat, ue)
379    USE wind_mod
380    REAL(rstd) :: dulon(iim*jjm,llm)
381    REAL(rstd) :: dulat(iim*jjm,llm)
382    REAL(rstd) :: ue(3*iim*jjm,llm)
383    REAL(rstd) :: duc(iim*jjm,llm,3)
384    REAL(rstd) :: dt2, due
385    INTEGER :: i,j,ij,l
386    ! Reconstruct wind tendencies at edges and add to normal wind
387    CALL compute_wind_centered_from_lonlat_compound(dulon,dulat,duc)
388    dt2=.5*physics_inout%dt_phys
389    DO l=ll_begin,ll_end
390      DO j=jj_begin,jj_end
391        DO i=ii_begin,ii_end
392          ij=(j-1)*iim+i
393          due = sum( (duc(ij,l,:) + duc(ij+t_right,l,:))*ep_e(ij+u_right,:) )
394          ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due
395
396          due = sum( (duc(ij,l,:) + duc(ij+t_lup,l,:))*ep_e(ij+u_lup,:) )
397          ue(ij+u_lup,l)=ue(ij+u_lup,l) + dt2*due
398
399          due = sum( (duc(ij,l,:) + duc(ij+t_ldown,l,:))*ep_e(ij+u_ldown,:) )
400          ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + dt2*due
401        ENDDO
402      ENDDO
403    ENDDO
404  END SUBROUTINE compute_update_velocity
405
406END MODULE physics_mod
Note: See TracBrowser for help on using the repository browser.