source: codes/icosagcm/devel/src/diagnostics/check_conserve.f90 @ 913

Last change on this file since 913 was 913, checked in by dubos, 5 years ago

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