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

Last change on this file since 171 was 160, checked in by dubos, 11 years ago

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