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

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