Changeset 714 for codes/icosagcm/devel/src/diagnostics/check_conserve.f90

Timestamp:

08/03/18 16:53:37 (6 years ago)

Author:

dubos

Message:

devel : backported from trunk commits r607,r648,r649,r667,r668,r669,r706

File:

• codes/icosagcm/devel/src/diagnostics/check_conserve.f90 (modified) (16 diffs)

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

@@ -14 +14 @@
   
   REAL(rstd),SAVE :: AAM_mass, AAM_mass_old, AAM_vel, AAM_vel_old, &
+       AAM_velp, AAM_velp_old, AAM_velm, AAM_velm_old, &
        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,rmsvtot0
   !$OMP THREADPRIVATE(check_type, mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0)
@@ -66 +68 @@
     INTEGER :: ind,ierr
     REAL(rstd) :: mtot, angtot, rmsdpdt
-    REAL(rstd) :: etot, stot, ang_mass, ang_vel, rmsvtot, ztot
+    REAL(rstd) :: etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsvtot, ztot
 
     CALL transfert_request(f_ue,req_e1_vect)
@@ -84 +86 @@
     CALL check_PV(ztot) 
     CALL exner(f_ps,f_p,f_pks,f_pk)
-    CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, rmsvtot) 
+    CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsvtot) 
 
     IF (is_master) THEN ! is_master = is_omp_master && is_mpi_master
@@ -90 +92 @@
        AAM_mass = ang_mass
        AAM_vel = ang_vel
+       AAM_velp = ang_velp
+       AAM_velm = ang_velm
        angtot  = ang_mass + ang_vel
        IF ( it == itau0  ) THEN 
@@ -99 +103 @@
           AAM_mass_old=AAM_mass
           AAM_vel_old=AAM_vel
+          AAM_velp_old=AAM_velp
+          AAM_velm_old=AAM_velm
           AAM_mass_source=0.
           AAM_vel_source=0.
+          AAM_vel_plus_source=0.
+          AAM_vel_minus_source=0.
        END IF
        rmsvtot=SQRT(rmsvtot/mtot) 
@@ -127 +135 @@
           WRITE(*,'(A,6E12.4)') 'AAM_vel  : time,old,new,dissip,dyn,phys', dt*it, AAM_vel_old, AAM_vel, &
                AAM_vel_source(AAM_dissip), AAM_vel_source(AAM_dyn), AAM_vel_source(AAM_phys)
-          WRITE(*,'(A,6E12.4)') 'AAM_tot  : time,old,new,dissip,dyn,phys', dt*it, &
-               AAM_vel_old+AAM_mass_old, AAM_vel+AAM_mass, &
-               AAM_vel_source(AAM_dissip)+AAM_mass_source(AAM_dissip), &
-               AAM_vel_source(AAM_dyn)   +AAM_mass_source(AAM_dyn), &
-               AAM_vel_source(AAM_phys)  +AAM_mass_source(AAM_phys)
+          WRITE(*,'(A,6E12.4)') 'AAM_vel+ : time,old,new,dissip,dyn,phys', dt*it, AAM_velp_old, AAM_velp, &
+               AAM_vel_plus_source(AAM_dissip), AAM_vel_plus_source(AAM_dyn), AAM_vel_plus_source(AAM_phys)
+          WRITE(*,'(A,6E12.4)') 'AAM_vel- : time,old,new,dissip,dyn,phys', dt*it, AAM_velm_old, AAM_velm, &
+               AAM_vel_minus_source(AAM_dissip), AAM_vel_minus_source(AAM_dyn), AAM_vel_minus_source(AAM_phys)
+          WRITE(*,'(A,6E12.4)') 'AAM_dyn budget mass + vel:', AAM_mass_source(AAM_dyn) + AAM_vel_source(AAM_dyn)
           AAM_mass_old=AAM_mass
           AAM_vel_old=AAM_vel
+          AAM_velp_old=AAM_velp
+          AAM_velm_old=AAM_velm
           AAM_mass_source=0.
           AAM_vel_source=0.
+          AAM_vel_minus_source=0.
+          AAM_vel_plus_source=0.
        END IF
     END IF
@@ -160 +172 @@
     REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 
     INTEGER::ind,ierr
-    REAL(rstd) :: mtot, ztot, rmsdpdt, etot,stot,rmsv, ang_mass, ang_vel
+    REAL(rstd) :: mtot, ztot, rmsdpdt, etot,stot,rmsv, ang_mass, ang_vel, ang_velp, ang_velm
     
     IF(check_type == check_detailed) THEN
@@ -175 +187 @@
        END DO
        CALL exner(f_ps,f_p,f_pks,f_pk)
-       CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, rmsv) 
+       CALL check_energy(f_ue,f_theta_rhodz,f_phis, etot, stot, ang_mass, ang_vel, ang_velp, ang_velm, rmsv) 
        
        IF (is_master) THEN
@@ -186 +198 @@
           AAM_mass_source(tag) = AAM_mass_source(tag) + ang_mass - AAM_mass
           AAM_vel_source(tag) = AAM_vel_source(tag) + ang_vel - AAM_vel
+          AAM_vel_plus_source(tag) = AAM_vel_plus_source(tag) + ang_velp - AAM_velp
+          AAM_vel_minus_source(tag) = AAM_vel_minus_source(tag) + ang_velm - AAM_velm
+          !if ((ang_vel - AAM_vel) .gt. 0.) then
+          !  AAM_vel_plus_source(tag) = AAM_vel_plus_source(tag) + ang_vel - AAM_vel
+          !else
+          !  AAM_vel_minus_source(tag) = AAM_vel_minus_source(tag) + ang_vel - AAM_vel
+          !endif
           AAM_mass = ang_mass
           AAM_vel = ang_vel
+          AAM_velp = ang_velp
+          AAM_velm = ang_velm
        END IF
     END IF
@@ -252 +273 @@
 
   SUBROUTINE check_energy(f_ue,f_theta_rhodz,f_phis, etot, &
-                      stot, AAM_mass_tot, AAM_vel_tot, rmsvtot)  
+                      stot, AAM_mass_tot, AAM_vel_tot, AAM_velp_tot, AAM_velm_tot, rmsvtot)  
   USE pression_mod
   USE vorticity_mod
@@ -258 +279 @@
     TYPE(t_field), POINTER :: f_theta_rhodz(:)
     TYPE(t_field), POINTER :: f_phis(:)
-    REAL(rstd),INTENT(OUT) :: etot, stot, AAM_mass_tot, AAM_vel_tot, rmsvtot
+    REAL(rstd),INTENT(OUT) :: etot, stot, AAM_mass_tot, AAM_vel_tot, AAM_velp_tot, AAM_velm_tot, rmsvtot
     REAL(rstd), POINTER :: ue(:,:)
     REAL(rstd), POINTER :: p(:,:) 
@@ -265 +286 @@
     REAL(rstd), POINTER :: phis(:) 
     REAL(rstd), POINTER :: rhodz(:,:) 
-    REAL(rstd) :: e, s, AAM_mass, AAM_vel, rmsv
+    REAL(rstd) :: e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv
     INTEGER :: ind
 
     e = 0 ; s = 0 ; AAM_mass = 0 ; AAM_vel=0 ; rmsv = 0
+    AAM_velp = 0 ; AAM_velm = 0
 
     DO ind=1,ndomain 
@@ -281 +303 @@
       phis=f_phis(ind) 
       CALL compute_energy(ind,ue,p,rhodz,theta_rhodz(:,:,1),pk,phis, &
-           e, s, AAM_mass, AAM_vel, rmsv)
+           e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv)
     END DO 
 
@@ -289 +311 @@
     CALL global_sum(AAM_mass, AAM_mass_tot)
     CALL global_sum(AAM_vel, AAM_vel_tot)
+    CALL global_sum(AAM_velp, AAM_velp_tot)
+    CALL global_sum(AAM_velm, AAM_velm_tot)
   END SUBROUTINE check_energy
 
-  SUBROUTINE compute_energy(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, AAM_mass, AAM_vel, rmsv)
+  SUBROUTINE compute_energy(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv)
     USE disvert_mod
     USE wind_mod
@@ -302 +326 @@
     REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
     REAL(rstd),INTENT(IN) :: phis(iim*jjm)
-    REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, rmsv
+    REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, AAM_velp, AAM_velm, rmsv
     
     REAL(rstd) :: ucenter(iim*jjm,llm,3)
@@ -335 +359 @@
                 AAM_mass = AAM_mass + rad*cos(lat)*mass_ij*radomeg*cos(lat)
                 AAM_vel  = AAM_vel  + rad*cos(lat)*mass_ij*ulon(ij,l)
+                if (ulon(ij,l) > 0.) then
+                  AAM_velp  = AAM_velp  + rad*cos(lat)*mass_ij*ulon(ij,l)
+                else
+                  AAM_velm  = AAM_velm  + rad*cos(lat)*mass_ij*ulon(ij,l)
+                endif
              END IF
           END DO