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

Last change on this file since 176 was 172, checked in by ymipsl, 11 years ago

doesn't work in sequential mode
=> Add MPI_REDUCE wrapper

YM
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    IF (using_mpi) THEN
129      CALL MPI_REDUCE(mloc, mtot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
130      CALL MPI_REDUCE(rmsloc, rmsdpdt, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
131    ENDIF
132   
133  END SUBROUTINE check_mass_conserve
134
135!---------------------------------------------------------------------
136
137  SUBROUTINE check_en(f_ue,f_theta_rhodz,f_phis) 
138  USE icosa
139  USE pression_mod
140  USE vorticity_mod
141  IMPLICIT NONE
142    TYPE(t_field), POINTER :: f_ue(:)
143    TYPE(t_field), POINTER :: f_theta_rhodz(:)
144    TYPE(t_field), POINTER :: f_phis(:)
145    REAL(rstd), POINTER :: ue(:,:)
146    REAL(rstd), POINTER :: p(:,:) 
147    REAL(rstd), POINTER :: theta_rhodz(:,:) 
148    REAL(rstd), POINTER :: pk(:,:) 
149    REAL(rstd), POINTER :: phis(:) 
150    REAL(rstd), POINTER :: rhodz(:,:) 
151    INTEGER :: ind
152
153
154    DO ind=1,ndomain 
155      CALL swap_dimensions(ind)
156      CALL swap_geometry(ind)
157      ue=f_ue(ind) 
158      p=f_p(ind)
159      rhodz=f_rhodz(ind) 
160      theta_rhodz=f_theta_rhodz(ind) 
161      pk=f_pk(ind) 
162      phis=f_phis(ind) 
163      CALL compute_en(ind,ue,p,rhodz,theta_rhodz,pk,phis)
164    END DO
165 
166  END SUBROUTINE check_en 
167
168  SUBROUTINE compute_en(ind,u,p,rhodz,theta_rhodz,pk,phis)
169  USE icosa
170  USE disvert_mod
171  USE wind_mod
172  IMPLICIT NONE
173    INTEGER,INTENT(IN)::ind 
174    REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm)
175    REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
176    REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
177    REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
178    REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
179    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
180    REAL(rstd):: theta(iim*jjm,llm)
181    REAL(rstd)::KE(iim*jjm,llm) 
182    REAL(rstd) :: ucenter(iim*jjm,3,llm)
183    REAL(rstd) :: ulon(iim*jjm,llm)
184    REAL(rstd) :: ulat(iim*jjm,llm)
185
186    REAL(rstd)::masse(iim*jjm,llm) 
187    REAL(rstd)::rad,radomeg,lat,lon
188    INTEGER :: i,j,ij,l,ij2
189
190 
191    DO l = 1, llm
192      DO j=jj_begin,jj_end
193        DO i=ii_begin,ii_end
194          ij=(j-1)*iim+i
195          masse(ij,l) = rhodz(ij,l)*Ai(ij) 
196          theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
197          IF (domain(ind)%own(i,j)) THEN
198            stot = stot + Ai(ij)*theta_rhodz(ij,l)
199          END IF
200        END DO
201      END DO
202    END DO
203
204    DO l = 1, llm
205      DO j=jj_begin,jj_end
206        DO i=ii_begin,ii_end
207          ij=(j-1)*iim+i
208          IF (domain(ind)%own(i,j)) THEN
209            KE(ij,l)=   1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +   &
210                                      le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +           &
211                                      le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +           &
212                                      le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +        &
213                                      le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +     &
214                                      le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
215            rmsv = rmsv + 2*masse(ij,l)*KE(ij,l) 
216          END IF
217        END DO
218      END DO
219    END DO 
220
221    DO l = 1, llm
222      DO j=jj_begin,jj_end
223        DO i=ii_begin,ii_end
224          ij=(j-1)*iim+i
225          IF (domain(ind)%own(i,j)) THEN
226            etot = etot + masse(ij,l)*(phis(ij)+theta(ij,l)*pk(ij,l)+KE(ij,l))
227          END IF
228        END DO
229      END DO
230    END DO 
231
232!------------------------------ ANGULAR VELOCITY
233    CALL compute_wind_centered(u,ucenter)
234    CALL compute_wind_centered_lonlat_compound(ucenter,ulon,ulat)
235    rad = radius 
236    radomeg = rad*omega
237
238    DO l = 1, llm
239      DO j=jj_begin,jj_end
240        DO i=ii_begin,ii_end
241          ij=(j-1)*iim+i
242          IF (domain(ind)%own(i,j)) THEN
243            CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 
244            ang=ang+rad*cos(lat)*masse(ij,l)*(ulon(ij,l)+ radomeg*cos(lat))
245          END IF
246        END DO
247      END DO
248    END DO   
249         
250  END SUBROUTINE compute_en
251
252!---------------------------------------------------------------------
253
254  SUBROUTINE check_PV 
255  USE icosa
256  IMPLICIT NONE
257     REAL(rstd), POINTER :: vort(:,:)
258     REAL(rstd), POINTER :: rhodz(:,:) 
259     INTEGER :: ind
260
261     DO ind=1,ndomain 
262       CALL swap_dimensions(ind)
263       CALL swap_geometry(ind)
264       vort=f_vort(ind)
265       rhodz=f_rhodz(ind) 
266       CALL compute_PV(vort,rhodz) 
267     ENDDO
268 
269  END SUBROUTINE check_PV 
270
271  SUBROUTINE compute_PV(vort,rhodz) 
272  USE icosa
273  USE disvert_mod
274  IMPLICIT NONE
275    REAL(rstd),INTENT(IN) :: vort(2*iim*jjm,llm)
276    REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 
277    REAL(rstd)::qv1,qv2 
278    REAL(rstd)::hv1,hv2 
279    INTEGER :: i,j,ij,l,ij2
280
281    hv1 = 0.0 ; hv2 = 0.0 
282   
283    DO l = 1,llm
284      DO j=jj_begin+1,jj_end-1
285        DO i=ii_begin+1,ii_end-1
286          ij=(j-1)*iim+i
287          hv1 =  Riv2(ij,vup)        * rhodz(ij,l)           &
288            + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
289            + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
290   
291          qv1 = ( vort(ij+z_up,l)+fv(ij+z_up) )/hv1
292         
293          hv2  =  Riv2(ij,vdown)     * rhodz(ij,l)          &
294            + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
295            + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
296                     
297          qv2 =( vort(ij+z_down,l)+fv(ij+z_down) )/hv2
298
299          ztot = ztot + (qv1*qv1*hv1 + qv2*qv2*hv2) 
300
301        ENDDO
302      ENDDO
303    ENDDO
304
305  END SUBROUTINE compute_PV
306!---------------------------------------------------------------------
307
308   SUBROUTINE compute_rhodz(p,rhodz) 
309   USE icosa 
310   IMPLICIT NONE
311     REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
312     REAL(rstd),INTENT(OUT):: rhodz(iim*jjm,llm)
313     REAL(rstd) ::mass(iim*jjm,llm)
314     INTEGER :: i,j,ij,l 
315       
316     DO l = 1, llm
317      DO j=jj_begin,jj_end
318       DO i=ii_begin,ii_end
319         ij=(j-1)*iim+i
320         mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g
321         rhodz(ij,l) = mass(ij,l) / Ai(ij)
322       ENDDO
323      ENDDO
324     ENDDO
325
326   END SUBROUTINE compute_rhodz
327!====================================================================
328
329
330 END MODULE check_conserve_mod
Note: See TracBrowser for help on using the repository browser.