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

Last change on this file since 609 was 609, checked in by dubos, 6 years ago

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