Changeset 902 for codes/icosagcm/trunk/src

Timestamp:

06/13/19 16:45:43 (5 years ago)

Author:

adurocher

Message:

trunk : Check tracer mass conservation in check_conserve

Location:

codes/icosagcm/trunk/src

Files:

• diagnostics/check_conserve.f90 (modified) (7 diffs)
• time/timeloop_gcm.f90 (modified) (3 diffs)

codes/icosagcm/trunk/src/diagnostics/check_conserve.f90

@@ -17 +17 @@
        AAM_mass_source(3), AAM_vel_source(3) ! read/written only IF is_master
   REAL(rstd),SAVE :: AAM_vel_plus_source(3), AAM_vel_minus_source(3)
-  REAL(rstd),SAVE :: mtot0,ztot0,etot0,angtot0,stot0
+  REAL(rstd),SAVE :: mtot0,mqtot0,ztot0,etot0,angtot0,stot0
   !$OMP THREADPRIVATE(check_type, mtot0,ztot0,etot0,angtot0,stot0)
   
@@ -51 +51 @@
 !--------------------------------- Basic check --------------------------------
 
-  SUBROUTINE check_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it) 
+  SUBROUTINE check_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,f_q,it) 
   USE pression_mod
   USE vorticity_mod
@@ -63 +63 @@
     TYPE(t_field),POINTER :: f_theta_rhodz(:)
     TYPE(t_field),POINTER :: f_phis(:)
+    TYPE(t_field),POINTER :: f_q(:)
     INTEGER, INTENT(IN) :: it
 
     REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 
     INTEGER :: ind
-    REAL(rstd) :: mtot, angtot, rmsdpdt
+    REAL(rstd) :: mtot, mqtot, angtot, rmsdpdt
     REAL(rstd) :: etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsvtot, ztot
 
@@ -83 +84 @@
 
     CALL vorticity(f_ue,f_vort)
+    CALL check_qmass_conserve(f_q,f_rhodz,mqtot)
     CALL check_mass_conserve(f_ps,f_dps,mtot,rmsdpdt)
     CALL check_PV(ztot) 
@@ -98 +100 @@
           ztot0 = ztot
           mtot0 = mtot
+          mqtot0 = mqtot
           etot0 = etot 
           angtot0  = angtot
@@ -111 +114 @@
        END IF
        rmsvtot=SQRT(rmsvtot/mtot) 
-       ztot=ztot/ztot0-1. ; mtot=mtot/mtot0-1. 
+       ztot=ztot/ztot0-1. ; mtot=mtot/mtot0-1.
+       mqtot=(mqtot+1E-100)/(mqtot0+1e-100)-1.
        etot=etot/etot0-1. ; angtot=angtot/angtot0-1. ; stot=stot/stot0-1. 
        rmsdpdt= daysec*1.e-2*sqrt(rmsdpdt/ncell_glo)  
        
        OPEN(134,file="checkconsicosa.txt",position='append')
-       WRITE(134,4000)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
-       WRITE(134,*)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot 
+       WRITE(134,4000)mtot,mqtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
+       WRITE(134,*)mtot,mqtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot 
        WRITE(134,*)"=================================================="
-       WRITE(*,4000)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
+       WRITE(*,4000)mtot,mqtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
        
-4000   FORMAT(10x,'masse',5x,'rmsdpdt',5x,'energie',5x,'enstrophie'  &
+4000   FORMAT(10x,'masse',5x,'advec mass',5x,'rmsdpdt',6x,'energie',3x,'enstrophie'  &
             ,5x,'entropie',5x,'rmsv',5x,'mt.ang',/,'GLOB  '          &
-            ,e10.3,e13.6,5e13.3/)     
+            ,2e10.3,e13.5,5e13.3/)     
        CLOSE(134)
 
@@ -269 +273 @@
   END SUBROUTINE check_mass_conserve
 
+  SUBROUTINE check_qmass_conserve(f_q,f_rhodz,mqtot)
+  USE mpi_mod
+  USE mpipara
+  USE transfert_omp_mod
+  USE omp_para
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_q(:)
+    TYPE(t_field),POINTER :: f_rhodz(:)
+    REAL(rstd),POINTER :: q(:,:,:) 
+    REAL(rstd),POINTER :: rhodz(:,:) 
+    REAL(rstd), INTENT(OUT) :: mqtot
+
+    INTEGER :: ind,i,j,ij,l 
+    REAL :: mloc
+
+    mloc=0.0
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      q=f_q(ind)
+      rhodz=f_rhodz(ind)
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          IF (domain(ind)%own(i,j)) THEN
+            DO l=1, llm
+             mloc=mloc + sum( q(ij,l,:)*Ai(ij)*rhodz(ij,l) )
+            END DO
+          ENDIF
+        ENDDO
+      ENDDO
+    ENDDO
+    
+    CALL global_sum(mloc, mqtot)
+  END SUBROUTINE check_qmass_conserve
+
 !---------------------------------------------------------------------
 

codes/icosagcm/trunk/src/time/timeloop_gcm.f90

@@ -228 +228 @@
     END IF
 
-    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,itau0)  
+    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,f_q,itau0)  
 
     Call trace_on
@@ -362 +362 @@
           CALL check_conserve_detailed(it, AAM_dyn, &
                f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
-          CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it) 
+          CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,f_q,it) 
        ENDIF
 
@@ -379 +379 @@
     CALL check_conserve_detailed(it, AAM_dyn, &
          f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
-    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)  
+    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,f_q,it)  
     
     !$OMP MASTER