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

Last change on this file since 156 was 151, checked in by ymipsl, 11 years ago

Implementation of mixte parallelism MPI/OpenMP into src directory

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