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

Last change on this file since 364 was 326, checked in by dubos, 9 years ago

Detailed energy/AAM diagnostics
File size: 12.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    INTEGER, PARAMETER :: check_basic=1, check_detailed=2
11    INTEGER, SAVE :: check_type
12  PUBLIC init_check_conserve, check_conserve, check_conserve_detailed
13    REAL(rstd),SAVE:: mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0
14!$OMP THREADPRIVATE(check_type, mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0)
15     
16CONTAINS 
17
18!---------------------------------------------------------------------
19
20  SUBROUTINE init_check_conserve
21    USE icosa
22    USE getin_mod
23    USE omp_para, ONLY : is_master
24    IMPLICIT NONE
25    CHARACTER(LEN=255) :: check_type_str
26    CALL allocate_field(f_pk,field_t,type_real,llm)
27    CALL allocate_field(f_p,field_t,type_real,llm+1)
28    CALL allocate_field(f_pks,field_t,type_real)
29    CALL allocate_field(f_rhodz,field_t,type_real,llm)
30    CALL allocate_field(f_vort,field_z,type_real,llm)
31   
32    check_type_str='basic'
33    CALL getin("check_conservation",check_type_str)
34    SELECT CASE(TRIM(check_type_str))
35    CASE('basic')
36       check_type=check_basic
37    CASE('detailed')
38       check_type=check_detailed
39    CASE DEFAULT
40       IF(is_master) PRINT*,'Bad selector <', TRIM(check_type_str), &
41            '> for variable check_conservation : options are <basic>, <detailed>'
42       STOP
43    END SELECT
44  END SUBROUTINE init_check_conserve
45
46!--------------------------------- Basic check --------------------------------
47
48  SUBROUTINE check_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it) 
49  USE icosa 
50  USE pression_mod
51  USE vorticity_mod
52  USE caldyn_gcm_mod
53  USE exner_mod 
54  USE mpipara, ONLY : is_mpi_master, comm_icosa
55  USE omp_para
56  IMPLICIT NONE
57    TYPE(t_field),POINTER :: f_ps(:)
58    TYPE(t_field),POINTER :: f_dps(:)
59    TYPE(t_field),POINTER :: f_ue(:)
60    TYPE(t_field),POINTER :: f_theta_rhodz(:)
61    TYPE(t_field),POINTER :: f_phis(:)
62    INTEGER, INTENT(IN) :: it
63
64    REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 
65    INTEGER::ind,ierr
66    REAL(rstd) :: mtot, angtot, rmsdpdt
67    REAL(rstd) :: etot, stot, ang_mass, ang_vit, rmsvtot, ztot
68
69    CALL transfert_request(f_ue,req_e1_vect)
70    CALL pression(f_ps,f_p)
71
72    DO ind=1,ndomain 
73      IF (.NOT. assigned_domain(ind)) CYCLE
74      CALL swap_dimensions(ind)
75      CALL swap_geometry(ind)
76      p=f_p(ind) 
77      rhodz=f_rhodz(ind)
78      CALL compute_rhodz(p,rhodz) 
79    END DO
80
81    CALL vorticity(f_ue,f_vort)
82    CALL check_mass_conserve(f_ps,f_dps,mtot,rmsdpdt)
83    CALL check_PV(ztot) 
84    CALL exner(f_ps,f_p,f_pks,f_pk)
85    CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vit, rmsvtot) 
86    angtot  = ang_mass + ang_vit 
87
88    IF (is_master) THEN
89       IF ( it == itau0  ) THEN
90          ztot0 = ztot
91          mtot0 = mtot
92          etot0 = etot 
93          angtot0  = angtot
94          stot0 = stot 
95       END IF
96
97       rmsvtot=SQRT(rmsvtot/mtot) 
98       ztot=ztot/ztot0-1. ; mtot=mtot/mtot0-1. 
99       etot=etot/etot0-1. ; angtot=angtot/angtot0-1. ; stot=stot/stot0-1. 
100       rmsdpdt= daysec*1.e-2*sqrt(rmsdpdt/ncell_glo) 
101       
102       OPEN(134,file="checkconsicosa.txt",position='append')
103       WRITE(134,4000)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
104       WRITE(134,*)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot 
105       WRITE(134,*)"=================================================="
106       WRITE(*,4000)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
107       
1084000   FORMAT(10x,'masse',5x,'rmsdpdt',5x,'energie',5x,'enstrophie'  &
109            ,5x,'entropie',5x,'rmsv',5x,'mt.ang',/,'GLOB  '                &
110            ,e10.3,e13.6,5e13.3/)     
111       close(134)
112     
113    END IF
114  END SUBROUTINE check_conserve
115
116!------------------------ Detailed check (Lauritzen et al., 2014) -----------------------------
117
118  SUBROUTINE check_conserve_detailed(tag, f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it) 
119    USE icosa 
120    USE pression_mod
121    USE vorticity_mod
122    USE caldyn_gcm_mod
123    USE exner_mod 
124    USE mpipara, ONLY : is_mpi_root, comm_icosa
125    IMPLICIT NONE
126    CHARACTER(LEN=*), INTENT(IN) :: tag
127    TYPE(t_field),POINTER :: f_ps(:)
128    TYPE(t_field),POINTER :: f_dps(:)
129    TYPE(t_field),POINTER :: f_ue(:)
130    TYPE(t_field),POINTER :: f_theta_rhodz(:)
131    TYPE(t_field),POINTER :: f_phis(:)
132    INTEGER, INTENT(IN) ::it
133
134    REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 
135    INTEGER::ind,ierr
136    REAL(rstd) :: mtot, ztot, rmsdpdt, etot,stot,rmsv, AAM_mass, AAM_vel
137   
138    IF(check_type == check_detailed) THEN
139       etot=0.0; stot=0.0; rmsv=0.0 
140       AAM_vel=0. ; AAM_vel=0.
141       ztot = 0.0 
142
143       CALL transfert_request(f_ue,req_e1_vect)
144       CALL pression(f_ps,f_p)
145       DO ind=1,ndomain 
146          IF (.NOT. assigned_domain(ind)) CYCLE
147          CALL swap_dimensions(ind)
148          CALL swap_geometry(ind)
149          p=f_p(ind) 
150          rhodz=f_rhodz(ind)
151          CALL compute_rhodz(p,rhodz) 
152       END DO
153       CALL exner(f_ps,f_p,f_pks,f_pk)
154       CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, AAM_mass, AAM_vel, rmsv) 
155       
156       IF (is_mpi_root) THEN 
157          !$OMP MASTER       
158          IF ( it == itau0  ) Then
159             ztot0 = ztot
160             mtot0 = mtot
161             etot0 = etot 
162             angtot0  = AAM_mass + AAM_vel
163             stot0 = stot 
164          END IF
165         
166          PRINT *, TRIM(tag), AAM_mass, AAM_vel
167          !$OMP END MASTER
168       END IF
169    END IF
170  END SUBROUTINE check_conserve_detailed
171 
172!---------------------------------------------------------------------
173
174  SUBROUTINE global_sum(loc_sum, glob_sum)
175    USE mpi_mod
176    USE mpipara
177    USE transfert_omp_mod
178    IMPLICIT NONE
179    REAL(rstd), INTENT(IN) :: loc_sum
180    REAL(rstd), INTENT(OUT) :: glob_sum
181    REAL(rstd) :: sum
182
183    CALL reduce_sum_omp(loc_sum, sum)
184    IF (using_mpi) THEN
185!$OMP MASTER   
186      CALL MPI_REDUCE(sum, glob_sum, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
187!$OMP END MASTER
188    ELSE
189      glob_sum=loc_sum
190    ENDIF
191  END SUBROUTINE global_sum
192
193  SUBROUTINE check_mass_conserve(f_ps,f_dps,mtot,rmsdpdt)
194  USE mpi_mod
195  USE mpipara
196  USE transfert_omp_mod
197  USE icosa
198  USE omp_para
199  IMPLICIT NONE
200    TYPE(t_field),POINTER :: f_ps(:)
201    TYPE(t_field),POINTER :: f_dps(:)
202    REAL(rstd),POINTER :: ps(:),dps(:) 
203    REAL(rstd), INTENT(OUT) :: mtot, rmsdpdt
204
205    INTEGER :: ind,i,j,ij 
206    REAL :: mloc, rmsloc
207    REAL :: mloc_mpi, rmsloc_mpi
208
209    mloc=0.0; rmsloc=0.0
210    DO ind=1,ndomain
211      IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
212      CALL swap_dimensions(ind)
213      CALL swap_geometry(ind)
214      dps=f_dps(ind) 
215      ps=f_ps(ind)
216      DO j=jj_begin,jj_end
217        DO i=ii_begin,ii_end
218          ij=(j-1)*iim+i
219          IF (domain(ind)%own(i,j)) THEN
220             mloc=mloc+ps(ij)*Ai(ij)
221             rmsloc=rmsloc+dps(ij)*dps(ij)
222          ENDIF
223        ENDDO
224      ENDDO
225    ENDDO
226   
227    CALL global_sum(mloc, mtot)
228    CALL global_sum(rmsloc, rmsdpdt)
229  END SUBROUTINE check_mass_conserve
230
231!---------------------------------------------------------------------
232
233  SUBROUTINE check_energy(f_ue,f_theta_rhodz,f_phis, etot, &
234                      stot, AAM_mass_tot, AAM_vel_tot, rmsvtot) 
235  USE icosa
236  USE pression_mod
237  USE vorticity_mod
238  IMPLICIT NONE
239    TYPE(t_field), POINTER :: f_ue(:)
240    TYPE(t_field), POINTER :: f_theta_rhodz(:)
241    TYPE(t_field), POINTER :: f_phis(:)
242    REAL(rstd),INTENT(OUT) :: etot, stot, AAM_mass_tot, AAM_vel_tot, rmsvtot
243    REAL(rstd), POINTER :: ue(:,:)
244    REAL(rstd), POINTER :: p(:,:) 
245    REAL(rstd), POINTER :: theta_rhodz(:,:) 
246    REAL(rstd), POINTER :: pk(:,:) 
247    REAL(rstd), POINTER :: phis(:) 
248    REAL(rstd), POINTER :: rhodz(:,:) 
249    REAL(rstd) :: e, s, AAM_mass, AAM_vel, rmsv
250    INTEGER :: ind
251
252    e = 0 ; s = 0 ; AAM_mass = 0 ; AAM_vel=0 ; rmsv = 0
253
254    DO ind=1,ndomain 
255      IF (.NOT. assigned_domain(ind)) CYCLE
256      CALL swap_dimensions(ind)
257      CALL swap_geometry(ind)
258      ue=f_ue(ind) 
259      p=f_p(ind)
260      rhodz=f_rhodz(ind) 
261      theta_rhodz=f_theta_rhodz(ind) 
262      pk=f_pk(ind) 
263      phis=f_phis(ind) 
264      CALL compute_energy(ind,ue,p,rhodz,theta_rhodz,pk,phis, &
265           e, s, AAM_mass, AAM_vel, rmsv)
266    END DO
267
268    CALL global_sum(e, etot)
269    CALL global_sum(s, stot)
270    CALL global_sum(rmsv, rmsvtot)
271    CALL global_sum(AAM_mass, AAM_mass_tot)
272    CALL global_sum(AAM_vel, AAM_vel_tot)
273  END SUBROUTINE check_energy
274
275  SUBROUTINE compute_energy(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, AAM_mass, AAM_vel, rmsv)
276    USE icosa
277    USE disvert_mod
278    USE wind_mod
279    USE omp_para
280    IMPLICIT NONE
281    INTEGER,INTENT(IN)::ind 
282    REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm)
283    REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
284    REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
285    REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
286    REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
287    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
288    REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, rmsv
289   
290    REAL(rstd) :: ucenter(iim*jjm,3,llm)
291    REAL(rstd) :: ulon(iim*jjm,llm)
292    REAL(rstd) :: ulat(iim*jjm,llm)
293   
294    REAL(rstd) :: mass_ij, theta_ij, KE_ij, rad,radomeg,lat
295    INTEGER :: i,j,ij,l
296   
297    CALL compute_wind_centered(u,ucenter)
298    CALL compute_wind_centered_lonlat_compound(ucenter,ulon,ulat)
299    rad = radius 
300    radomeg = rad*omega
301   
302    DO l = ll_begin,ll_end
303       DO j=jj_begin,jj_end
304          DO i=ii_begin,ii_end
305             ij=(j-1)*iim+i
306             IF (domain(ind)%own(i,j)) THEN
307                lat = lat_i(ij) 
308                mass_ij = rhodz(ij,l)*Ai(ij) 
309                theta_ij = theta_rhodz(ij,l)/rhodz(ij,l)
310                KE_ij = 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +   &
311                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +           &
312                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +           &
313                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +        &
314                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +     &
315                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
316                rmsv = rmsv + 2*mass_ij*KE_ij
317                s    = s + Ai(ij)*theta_rhodz(ij,l)
318                e    = e + mass_ij*(phis(ij)+theta_ij*pk(ij,l)+KE_ij)
319                AAM_mass = AAM_mass + rad*cos(lat)*mass_ij*radomeg*cos(lat)
320                AAM_vel  = AAM_vel  + rad*cos(lat)*mass_ij*ulon(ij,l)
321             END IF
322          END DO
323       END DO
324    END DO
325   
326END SUBROUTINE compute_energy
327
328!---------------------------------------------------------------------
329
330  SUBROUTINE check_PV(ztot)
331  USE icosa
332  USE mpi_mod
333  USE mpipara
334  USE transfert_omp_mod
335  IMPLICIT NONE
336     REAL(rstd),INTENT(OUT) :: ztot
337     REAL(rstd), POINTER :: vort(:,:)
338     REAL(rstd), POINTER :: rhodz(:,:) 
339     INTEGER :: ind
340     REAL(rstd) :: z,z_mpi
341
342     z=0
343     DO ind=1,ndomain 
344       IF (.NOT. assigned_domain(ind)) CYCLE
345       CALL swap_dimensions(ind)
346       CALL swap_geometry(ind)
347       vort=f_vort(ind)
348       rhodz=f_rhodz(ind) 
349       CALL compute_PV(vort,rhodz,z) 
350     ENDDO
351
352    IF (using_mpi) THEN
353      CALL reduce_sum_omp(z, z_mpi)
354!$OMP MASTER   
355      CALL MPI_REDUCE(z_mpi, ztot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
356!$OMP END MASTER
357    ELSE
358      ztot=z
359    ENDIF
360 
361  END SUBROUTINE check_PV 
362
363  SUBROUTINE compute_PV(vort,rhodz,z) 
364  USE icosa
365  USE disvert_mod
366  USE omp_para
367  IMPLICIT NONE
368    REAL(rstd),INTENT(IN) :: vort(2*iim*jjm,llm)
369    REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 
370    REAL(rstd),INTENT(INOUT) :: z
371    REAL(rstd)::qv1,qv2 
372    REAL(rstd)::hv1,hv2 
373    INTEGER :: i,j,ij,l,ij2
374
375    hv1 = 0.0 ; hv2 = 0.0 
376   
377    DO l = ll_begin,ll_end
378      DO j=jj_begin+1,jj_end-1
379        DO i=ii_begin+1,ii_end-1
380          ij=(j-1)*iim+i
381          hv1 =  Riv2(ij,vup)        * rhodz(ij,l)           &
382            + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
383            + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
384   
385          qv1 = ( vort(ij+z_up,l)+fv(ij+z_up) )/hv1
386         
387          hv2  =  Riv2(ij,vdown)     * rhodz(ij,l)          &
388            + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
389            + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
390                     
391          qv2 =( vort(ij+z_down,l)+fv(ij+z_down) )/hv2
392
393          z = z + (qv1*qv1*hv1 + qv2*qv2*hv2) 
394
395        ENDDO
396      ENDDO
397     
398    ENDDO
399
400  END SUBROUTINE compute_PV
401!---------------------------------------------------------------------
402
403   SUBROUTINE compute_rhodz(p,rhodz) 
404   USE icosa 
405   IMPLICIT NONE
406     REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
407     REAL(rstd),INTENT(OUT):: rhodz(iim*jjm,llm)
408     REAL(rstd) ::mass(iim*jjm,llm)
409     INTEGER :: i,j,ij,l 
410       
411     DO l = 1, llm
412      DO j=jj_begin,jj_end
413       DO i=ii_begin,ii_end
414         ij=(j-1)*iim+i
415         mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g
416         rhodz(ij,l) = mass(ij,l) / Ai(ij)
417       ENDDO
418      ENDDO
419     ENDDO
420
421   END SUBROUTINE compute_rhodz
422!====================================================================
423
424
425 END MODULE check_conserve_mod
Note: See TracBrowser for help on using the repository browser.