source: codes/icosagcm/trunk/src/check_conserve.f90 @ 324

Last change on this file since 324 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: 11.4 KB
Line 
1MODULE check_conserve_mod
2  USE icosa 
3  IMPLICIT NONE
4
5  PRIVATE
6    TYPE(t_field),POINTER,SAVE :: f_pk(:), f_pks(:), f_p(:)
7    TYPE(t_field),POINTER,SAVE :: f_vort(:)
8    TYPE(t_field),POINTER,SAVE :: f_rhodz(:) 
9   
10  PUBLIC init_check_conserve, check_conserve 
11    REAL(rstd),SAVE:: mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0
12!$OMP THREADPRIVATE(mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0)
13     
14CONTAINS 
15
16!---------------------------------------------------------------------
17
18  SUBROUTINE init_check_conserve
19  USE icosa
20  IMPLICIT NONE
21    CALL allocate_field(f_pk,field_t,type_real,llm)
22    CALL allocate_field(f_p,field_t,type_real,llm+1)
23    CALL allocate_field(f_pks,field_t,type_real)
24    CALL allocate_field(f_rhodz,field_t,type_real,llm)
25    CALL allocate_field(f_vort,field_z,type_real,llm)
26  END SUBROUTINE init_check_conserve
27
28  SUBROUTINE check_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it) 
29  USE icosa 
30  USE pression_mod
31  USE vorticity_mod
32  USE caldyn_gcm_mod
33  USE exner_mod 
34  USE mpipara, ONLY : is_mpi_master, comm_icosa
35  USE omp_para
36  IMPLICIT NONE
37    TYPE(t_field),POINTER :: f_ps(:)
38    TYPE(t_field),POINTER :: f_dps(:)
39    TYPE(t_field),POINTER :: f_ue(:)
40    TYPE(t_field),POINTER :: f_theta_rhodz(:)
41    TYPE(t_field),POINTER :: f_phis(:)
42    INTEGER::it
43
44    REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 
45    INTEGER::ind,ierr
46    REAL(rstd) :: mtot, rmsdpdt
47    REAL(rstd) :: etot, stot, angtot, rmsvtot, ztot
48
49    CALL transfert_request(f_ue,req_e1_vect)
50    CALL pression(f_ps,f_p)
51
52    DO ind=1,ndomain 
53      IF (.NOT. assigned_domain(ind)) CYCLE
54      CALL swap_dimensions(ind)
55      CALL swap_geometry(ind)
56      p=f_p(ind) 
57      rhodz=f_rhodz(ind)
58      CALL compute_rhodz(p,rhodz) 
59    END DO
60
61    CALL vorticity(f_ue,f_vort)
62    CALL check_mass_conserve(f_ps,f_dps,mtot,rmsdpdt)
63    CALL check_PV(ztot) 
64    CALL exner(f_ps,f_p,f_pks,f_pk)
65    CALL check_EN(f_ue,f_theta_rhodz,f_phis, etot, stot, angtot, rmsvtot) 
66
67    IF (is_master) THEN
68       IF ( it == itau0  ) THEN
69          ztot0 = ztot
70          mtot0 = mtot
71          etot0 = etot 
72          angtot0  = angtot 
73          stot0 = stot 
74       END IF
75
76       rmsvtot=SQRT(rmsvtot/mtot) 
77       ztot=ztot/ztot0-1. ; mtot=mtot/mtot0-1. 
78       etot=etot/etot0-1. ; angtot=angtot/angtot0-1. ; stot=stot/stot0-1. 
79       rmsdpdt= daysec*1.e-2*sqrt(rmsdpdt/ncell_glo) 
80       
81       OPEN(134,file="checkconsicosa.txt",position='append')
82       WRITE(134,4000)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
83       WRITE(134,*)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot 
84       WRITE(134,*)"=================================================="
85       WRITE(*,4000)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
86       
874000   FORMAT(10x,'masse',5x,'rmsdpdt',5x,'energie',5x,'enstrophie'  &
88            ,5x,'entropie',5x,'rmsv',5x,'mt.ang',/,'GLOB  '                &
89            ,e10.3,e13.6,5e13.3/)     
90       close(134)
91     
92    END IF
93  END SUBROUTINE check_conserve
94 
95!---------------------------------------------------------------------
96
97  SUBROUTINE check_mass_conserve(f_ps,f_dps,mtot,rmsdpdt)
98  USE mpi_mod
99  USE mpipara
100  USE transfert_omp_mod
101  USE icosa
102  USE omp_para
103  IMPLICIT NONE
104    TYPE(t_field),POINTER :: f_ps(:)
105    TYPE(t_field),POINTER :: f_dps(:)
106    REAL(rstd),POINTER :: ps(:),dps(:) 
107    REAL(rstd), INTENT(OUT) :: mtot, rmsdpdt
108
109    INTEGER :: ind,i,j,ij 
110    REAL :: mloc, rmsloc
111    REAL :: mloc_mpi, rmsloc_mpi
112
113    mloc=0.0; rmsloc=0.0
114    DO ind=1,ndomain
115      IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
116      CALL swap_dimensions(ind)
117      CALL swap_geometry(ind)
118      dps=f_dps(ind) 
119      ps=f_ps(ind)
120      DO j=jj_begin,jj_end
121        DO i=ii_begin,ii_end
122          ij=(j-1)*iim+i
123          IF (domain(ind)%own(i,j)) THEN
124             mloc=mloc+ps(ij)*Ai(ij)
125             rmsloc=rmsloc+dps(ij)*dps(ij)
126          ENDIF
127        ENDDO
128      ENDDO
129    ENDDO
130
131    IF (using_mpi) THEN
132        CALL reduce_sum_omp(mloc, mloc_mpi)
133        CALL reduce_sum_omp(rmsloc, rmsloc_mpi)
134!$OMP MASTER   
135      CALL MPI_REDUCE(mloc_mpi, mtot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
136      CALL MPI_REDUCE(rmsloc_mpi, rmsdpdt, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
137!$OMP END MASTER
138    ELSE
139      mtot=mloc
140      rmsdpdt=rmsloc
141    ENDIF
142   
143  END SUBROUTINE check_mass_conserve
144
145!---------------------------------------------------------------------
146
147  SUBROUTINE check_en(f_ue,f_theta_rhodz,f_phis, etot, stot, angtot, rmsvtot) 
148  USE icosa
149  USE pression_mod
150  USE vorticity_mod
151  USE mpi_mod
152  USE mpipara
153  USE transfert_omp_mod
154  IMPLICIT NONE
155    TYPE(t_field), POINTER :: f_ue(:)
156    TYPE(t_field), POINTER :: f_theta_rhodz(:)
157    TYPE(t_field), POINTER :: f_phis(:)
158    REAL(rstd),INTENT(OUT) :: etot, stot, angtot, rmsvtot
159    REAL(rstd), POINTER :: ue(:,:)
160    REAL(rstd), POINTER :: p(:,:) 
161    REAL(rstd), POINTER :: theta_rhodz(:,:) 
162    REAL(rstd), POINTER :: pk(:,:) 
163    REAL(rstd), POINTER :: phis(:) 
164    REAL(rstd), POINTER :: rhodz(:,:) 
165    REAL(rstd) :: e, s, ang, rmsv
166    REAL(rstd) :: e_mpi, s_mpi, ang_mpi, rmsv_mpi
167    INTEGER :: ind
168
169    e = 0 ; s = 0 ; ang = 0 ; rmsv = 0
170
171    DO ind=1,ndomain 
172      IF (.NOT. assigned_domain(ind)) CYCLE
173      CALL swap_dimensions(ind)
174      CALL swap_geometry(ind)
175      ue=f_ue(ind) 
176      p=f_p(ind)
177      rhodz=f_rhodz(ind) 
178      theta_rhodz=f_theta_rhodz(ind) 
179      pk=f_pk(ind) 
180      phis=f_phis(ind) 
181      CALL compute_en(ind,ue,p,rhodz,theta_rhodz,pk,phis, e, s, ang, rmsv)
182    END DO
183
184    IF (using_mpi) THEN
185      CALL reduce_sum_omp(e, e_mpi)
186      CALL reduce_sum_omp(s, s_mpi)
187      CALL reduce_sum_omp(ang, ang_mpi)
188      CALL reduce_sum_omp(rmsv, rmsv_mpi)
189!$OMP MASTER   
190      CALL MPI_REDUCE(e_mpi, etot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
191      CALL MPI_REDUCE(s_mpi, stot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
192      CALL MPI_REDUCE(ang_mpi, angtot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
193      CALL MPI_REDUCE(rmsv_mpi, rmsvtot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
194!$OMP END MASTER
195    ELSE
196      etot=e
197      stot=s
198      angtot=ang
199      rmsvtot=rmsv
200    ENDIF
201
202  END SUBROUTINE check_en 
203
204  SUBROUTINE compute_en(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, ang, rmsv)
205  USE icosa
206  USE disvert_mod
207  USE wind_mod
208  USE omp_para
209  IMPLICIT NONE
210    INTEGER,INTENT(IN)::ind 
211    REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm)
212    REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
213    REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
214    REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
215    REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
216    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
217    REAL(rstd),INTENT(INOUT) :: e
218    REAL(rstd),INTENT(INOUT) :: s
219    REAL(rstd),INTENT(INOUT) :: ang
220    REAL(rstd),INTENT(INOUT) :: rmsv
221   
222    REAL(rstd):: theta(iim*jjm,llm)
223    REAL(rstd)::KE(iim*jjm,llm) 
224    REAL(rstd) :: ucenter(iim*jjm,3,llm)
225    REAL(rstd) :: ulon(iim*jjm,llm)
226    REAL(rstd) :: ulat(iim*jjm,llm)
227
228    REAL(rstd)::masse(iim*jjm,llm) 
229    REAL(rstd)::rad,radomeg,lat,lon
230    REAL(rstd) ::xyz(3)
231    INTEGER :: i,j,ij,l,ij2
232
233 
234    DO l = ll_begin, ll_end
235      DO j=jj_begin,jj_end
236        DO i=ii_begin,ii_end
237          ij=(j-1)*iim+i
238          masse(ij,l) = rhodz(ij,l)*Ai(ij) 
239          theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
240          IF (domain(ind)%own(i,j)) THEN
241            s = s + Ai(ij)*theta_rhodz(ij,l)
242          END IF
243        END DO
244      END DO
245    END DO
246
247    DO l = ll_begin,ll_end
248      DO j=jj_begin,jj_end
249        DO i=ii_begin,ii_end
250          ij=(j-1)*iim+i
251          IF (domain(ind)%own(i,j)) THEN
252            KE(ij,l)=   1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +   &
253                                      le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +           &
254                                      le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +           &
255                                      le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +        &
256                                      le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +     &
257                                      le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
258            rmsv = rmsv + 2*masse(ij,l)*KE(ij,l) 
259          END IF
260        END DO
261      END DO
262    END DO 
263
264    DO l = ll_begin,ll_end
265      DO j=jj_begin,jj_end
266        DO i=ii_begin,ii_end
267          ij=(j-1)*iim+i
268          IF (domain(ind)%own(i,j)) THEN
269            e = e + masse(ij,l)*(phis(ij)+theta(ij,l)*pk(ij,l)+KE(ij,l))
270          END IF
271        END DO
272      END DO
273    END DO 
274
275!------------------------------ ANGULAR VELOCITY
276    CALL compute_wind_centered(u,ucenter)
277    CALL compute_wind_centered_lonlat_compound(ucenter,ulon,ulat)
278    rad = radius 
279    radomeg = rad*omega
280
281    DO l = ll_begin,ll_end
282      DO j=jj_begin,jj_end
283        DO i=ii_begin,ii_end
284          ij=(j-1)*iim+i
285          IF (domain(ind)%own(i,j)) THEN
286            lat=lat_i(ij) 
287            ang=ang+rad*cos(lat)*masse(ij,l)*(ulon(ij,l)+ radomeg*cos(lat))
288          END IF
289        END DO
290      END DO
291    END DO   
292         
293  END SUBROUTINE compute_en
294
295!---------------------------------------------------------------------
296
297  SUBROUTINE check_PV(ztot)
298  USE icosa
299  USE mpi_mod
300  USE mpipara
301  USE transfert_omp_mod
302  IMPLICIT NONE
303     REAL(rstd),INTENT(OUT) :: ztot
304     REAL(rstd), POINTER :: vort(:,:)
305     REAL(rstd), POINTER :: rhodz(:,:) 
306     INTEGER :: ind
307     REAL(rstd) :: z,z_mpi
308
309     z=0
310     DO ind=1,ndomain 
311       IF (.NOT. assigned_domain(ind)) CYCLE
312       CALL swap_dimensions(ind)
313       CALL swap_geometry(ind)
314       vort=f_vort(ind)
315       rhodz=f_rhodz(ind) 
316       CALL compute_PV(vort,rhodz,z) 
317     ENDDO
318
319    IF (using_mpi) THEN
320      CALL reduce_sum_omp(z, z_mpi)
321!$OMP MASTER   
322      CALL MPI_REDUCE(z_mpi, ztot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
323!$OMP END MASTER
324    ELSE
325      ztot=z
326    ENDIF
327 
328  END SUBROUTINE check_PV 
329
330  SUBROUTINE compute_PV(vort,rhodz,z) 
331  USE icosa
332  USE disvert_mod
333  USE omp_para
334  IMPLICIT NONE
335    REAL(rstd),INTENT(IN) :: vort(2*iim*jjm,llm)
336    REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 
337    REAL(rstd),INTENT(INOUT) :: z
338    REAL(rstd)::qv1,qv2 
339    REAL(rstd)::hv1,hv2 
340    INTEGER :: i,j,ij,l,ij2
341
342    hv1 = 0.0 ; hv2 = 0.0 
343   
344    DO l = ll_begin,ll_end
345      DO j=jj_begin+1,jj_end-1
346        DO i=ii_begin+1,ii_end-1
347          ij=(j-1)*iim+i
348          hv1 =  Riv2(ij,vup)        * rhodz(ij,l)           &
349            + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
350            + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
351   
352          qv1 = ( vort(ij+z_up,l)+fv(ij+z_up) )/hv1
353         
354          hv2  =  Riv2(ij,vdown)     * rhodz(ij,l)          &
355            + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
356            + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
357                     
358          qv2 =( vort(ij+z_down,l)+fv(ij+z_down) )/hv2
359
360          z = z + (qv1*qv1*hv1 + qv2*qv2*hv2) 
361
362        ENDDO
363      ENDDO
364     
365    ENDDO
366
367  END SUBROUTINE compute_PV
368!---------------------------------------------------------------------
369
370   SUBROUTINE compute_rhodz(p,rhodz) 
371   USE icosa 
372   IMPLICIT NONE
373     REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
374     REAL(rstd),INTENT(OUT):: rhodz(iim*jjm,llm)
375     REAL(rstd) ::mass(iim*jjm,llm)
376     INTEGER :: i,j,ij,l 
377       
378     DO l = 1, llm
379      DO j=jj_begin,jj_end
380       DO i=ii_begin,ii_end
381         ij=(j-1)*iim+i
382         mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g
383         rhodz(ij,l) = mass(ij,l) / Ai(ij)
384       ENDDO
385      ENDDO
386     ENDDO
387
388   END SUBROUTINE compute_rhodz
389!====================================================================
390
391
392 END MODULE check_conserve_mod
Note: See TracBrowser for help on using the repository browser.