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

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

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

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