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

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

For intel compiler, bad inlining when using multithreading and high optimization level, causing crash.
Call has been replaced (temporarily) by corresponding code.

YM
File size: 10.1 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    REAL(rstd) ::xyz(3)
203    INTEGER :: i,j,ij,l,ij2
204
205 
206    DO l = 1, llm
207      DO j=jj_begin,jj_end
208        DO i=ii_begin,ii_end
209          ij=(j-1)*iim+i
210          masse(ij,l) = rhodz(ij,l)*Ai(ij) 
211          theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
212          IF (domain(ind)%own(i,j)) THEN
213            stot = stot + Ai(ij)*theta_rhodz(ij,l)
214          END IF
215        END DO
216      END DO
217    END DO
218
219    DO l = 1, llm
220      DO j=jj_begin,jj_end
221        DO i=ii_begin,ii_end
222          ij=(j-1)*iim+i
223          IF (domain(ind)%own(i,j)) THEN
224            KE(ij,l)=   1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +   &
225                                      le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +           &
226                                      le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +           &
227                                      le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +        &
228                                      le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +     &
229                                      le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
230            rmsv = rmsv + 2*masse(ij,l)*KE(ij,l) 
231          END IF
232        END DO
233      END DO
234    END DO 
235
236    DO l = 1, llm
237      DO j=jj_begin,jj_end
238        DO i=ii_begin,ii_end
239          ij=(j-1)*iim+i
240          IF (domain(ind)%own(i,j)) THEN
241            etot = etot + masse(ij,l)*(phis(ij)+theta(ij,l)*pk(ij,l)+KE(ij,l))
242          END IF
243        END DO
244      END DO
245    END DO 
246
247!------------------------------ ANGULAR VELOCITY
248    CALL compute_wind_centered(u,ucenter)
249    CALL compute_wind_centered_lonlat_compound(ucenter,ulon,ulat)
250    rad = radius 
251    radomeg = rad*omega
252
253    DO l = 1, llm
254      DO j=jj_begin,jj_end
255        DO i=ii_begin,ii_end
256          ij=(j-1)*iim+i
257          IF (domain(ind)%own(i,j)) THEN
258!            CALL xyz2lonlat(xyz_i(ij,:),lon,lat)
259! very bad inlining with intel compiler (>14 ?)
260             xyz(:)=xyz_i(ij,:)
261             xyz(:)=xyz(:)/sqrt(sum(xyz(:)**2))
262             lat=asin(xyz(3))
263             lon=atan2(xyz(2),xyz(1))
264 
265            ang=ang+rad*cos(lat)*masse(ij,l)*(ulon(ij,l)+ radomeg*cos(lat))
266          END IF
267        END DO
268      END DO
269    END DO   
270         
271  END SUBROUTINE compute_en
272
273!---------------------------------------------------------------------
274
275  SUBROUTINE check_PV 
276  USE icosa
277  IMPLICIT NONE
278     REAL(rstd), POINTER :: vort(:,:)
279     REAL(rstd), POINTER :: rhodz(:,:) 
280     INTEGER :: ind
281
282     DO ind=1,ndomain 
283       IF (.NOT. assigned_domain(ind)) CYCLE
284       CALL swap_dimensions(ind)
285       CALL swap_geometry(ind)
286       vort=f_vort(ind)
287       rhodz=f_rhodz(ind) 
288       CALL compute_PV(vort,rhodz) 
289     ENDDO
290 
291  END SUBROUTINE check_PV 
292
293  SUBROUTINE compute_PV(vort,rhodz) 
294  USE icosa
295  USE disvert_mod
296  IMPLICIT NONE
297    REAL(rstd),INTENT(IN) :: vort(2*iim*jjm,llm)
298    REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 
299    REAL(rstd)::qv1,qv2 
300    REAL(rstd)::hv1,hv2 
301    INTEGER :: i,j,ij,l,ij2
302
303    hv1 = 0.0 ; hv2 = 0.0 
304   
305    DO l = 1,llm
306      DO j=jj_begin+1,jj_end-1
307        DO i=ii_begin+1,ii_end-1
308          ij=(j-1)*iim+i
309          hv1 =  Riv2(ij,vup)        * rhodz(ij,l)           &
310            + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
311            + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
312   
313          qv1 = ( vort(ij+z_up,l)+fv(ij+z_up) )/hv1
314         
315          hv2  =  Riv2(ij,vdown)     * rhodz(ij,l)          &
316            + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
317            + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
318                     
319          qv2 =( vort(ij+z_down,l)+fv(ij+z_down) )/hv2
320
321          ztot = ztot + (qv1*qv1*hv1 + qv2*qv2*hv2) 
322
323        ENDDO
324      ENDDO
325    ENDDO
326
327  END SUBROUTINE compute_PV
328!---------------------------------------------------------------------
329
330   SUBROUTINE compute_rhodz(p,rhodz) 
331   USE icosa 
332   IMPLICIT NONE
333     REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
334     REAL(rstd),INTENT(OUT):: rhodz(iim*jjm,llm)
335     REAL(rstd) ::mass(iim*jjm,llm)
336     INTEGER :: i,j,ij,l 
337       
338     DO l = 1, llm
339      DO j=jj_begin,jj_end
340       DO i=ii_begin,ii_end
341         ij=(j-1)*iim+i
342         mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g
343         rhodz(ij,l) = mass(ij,l) / Ai(ij)
344       ENDDO
345      ENDDO
346     ENDDO
347
348   END SUBROUTINE compute_rhodz
349!====================================================================
350
351
352 END MODULE check_conserve_mod
Note: See TracBrowser for help on using the repository browser.