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

Last change on this file since 286 was 286, checked in by dubos, 10 years ago

Partial etat0 cleanup (removed calls to xyz2lonlat)
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 == itau0  ) 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            lat=lat_i(ij) 
259            ang=ang+rad*cos(lat)*masse(ij,l)*(ulon(ij,l)+ radomeg*cos(lat))
260          END IF
261        END DO
262      END DO
263    END DO   
264         
265  END SUBROUTINE compute_en
266
267!---------------------------------------------------------------------
268
269  SUBROUTINE check_PV 
270  USE icosa
271  IMPLICIT NONE
272     REAL(rstd), POINTER :: vort(:,:)
273     REAL(rstd), POINTER :: rhodz(:,:) 
274     INTEGER :: ind
275
276     DO ind=1,ndomain 
277       IF (.NOT. assigned_domain(ind)) CYCLE
278       CALL swap_dimensions(ind)
279       CALL swap_geometry(ind)
280       vort=f_vort(ind)
281       rhodz=f_rhodz(ind) 
282       CALL compute_PV(vort,rhodz) 
283     ENDDO
284 
285  END SUBROUTINE check_PV 
286
287  SUBROUTINE compute_PV(vort,rhodz) 
288  USE icosa
289  USE disvert_mod
290  IMPLICIT NONE
291    REAL(rstd),INTENT(IN) :: vort(2*iim*jjm,llm)
292    REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 
293    REAL(rstd)::qv1,qv2 
294    REAL(rstd)::hv1,hv2 
295    INTEGER :: i,j,ij,l,ij2
296
297    hv1 = 0.0 ; hv2 = 0.0 
298   
299    DO l = 1,llm
300      DO j=jj_begin+1,jj_end-1
301        DO i=ii_begin+1,ii_end-1
302          ij=(j-1)*iim+i
303          hv1 =  Riv2(ij,vup)        * rhodz(ij,l)           &
304            + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
305            + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
306   
307          qv1 = ( vort(ij+z_up,l)+fv(ij+z_up) )/hv1
308         
309          hv2  =  Riv2(ij,vdown)     * rhodz(ij,l)          &
310            + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
311            + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
312                     
313          qv2 =( vort(ij+z_down,l)+fv(ij+z_down) )/hv2
314
315          ztot = ztot + (qv1*qv1*hv1 + qv2*qv2*hv2) 
316
317        ENDDO
318      ENDDO
319    ENDDO
320
321  END SUBROUTINE compute_PV
322!---------------------------------------------------------------------
323
324   SUBROUTINE compute_rhodz(p,rhodz) 
325   USE icosa 
326   IMPLICIT NONE
327     REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
328     REAL(rstd),INTENT(OUT):: rhodz(iim*jjm,llm)
329     REAL(rstd) ::mass(iim*jjm,llm)
330     INTEGER :: i,j,ij,l 
331       
332     DO l = 1, llm
333      DO j=jj_begin,jj_end
334       DO i=ii_begin,ii_end
335         ij=(j-1)*iim+i
336         mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g
337         rhodz(ij,l) = mass(ij,l) / Ai(ij)
338       ENDDO
339      ENDDO
340     ENDDO
341
342   END SUBROUTINE compute_rhodz
343!====================================================================
344
345
346 END MODULE check_conserve_mod
Note: See TracBrowser for help on using the repository browser.