Changeset 3393 for trunk/LMDZ.MARS/libf/phymars/nonoro_gwd_mix_mod.F90

Timestamp:

Jul 4, 2024, 4:14:10 PM (5 months ago)

Author:

jliu

Message:

Update of non-orographic gravity waves mixing scheme. 1)mixing in potential
temperature are added. 2)all mixing in q,u,theta now are implemented by AR-1
algorithm. Tests (MY29,MY32-35) runs show that this implementation has limited
impact to the temperature/tides fields.

File:

• trunk/LMDZ.MARS/libf/phymars/nonoro_gwd_mix_mod.F90 (modified) (20 diffs)

trunk/LMDZ.MARS/libf/phymars/nonoro_gwd_mix_mod.F90

@@ -5 +5 @@
 REAL,allocatable,save :: du_eddymix_gwd(:,:) ! Zonal wind tendency due to GWD 
 REAL,allocatable,save :: dv_eddymix_gwd(:,:) ! Meridional wind tendency due to GWD
+REAL,allocatable,save :: dh_eddymix_gwd(:,:) ! PT tendency due to GWD
+REAL,allocatable,save :: dq_eddymix_gwd(:,:,:) ! tracers tendency due to GWD
 REAL,allocatable,save :: de_eddymix_rto(:,:) ! Meridional wind tendency due to GWD
 REAL,allocatable,save :: df_eddymix_flx(:,:) ! Meridional wind tendency due to GWD
@@ -11 +13 @@
 LOGICAL, save :: calljliu_gwimix ! flag for using non-orographic GW-induced mixing
 
-!$OMP THREADPRIVATE(du_eddymix_gwd,dv_eddymix_gwd,de_eddymix_rto,df_eddymix_flx,calljliu_gwimix) 
+!$OMP THREADPRIVATE(du_eddymix_gwd,dv_eddymix_gwd,dh_eddymix_gwd,dq_eddymix_gwd,de_eddymix_rto,df_eddymix_flx,calljliu_gwimix) 
 !,east_gwstress,west_gwstress)
 
@@ -17 +19 @@
 
       SUBROUTINE NONORO_GWD_MIX(ngrid,nlayer,DTIME,nq,cpnew, rnew, pp,  &
-                  zmax_therm, pt, pu, pv, pq, pdt, pdu, pdv, pdq,       &
+                  zmax_therm, pt, pu, pv, pq, pht, pdt, pdu, pdv, pdq, pdht, &
                   d_pq, d_t, d_u, d_v)
 
@@ -29 +31 @@
     !---------------------------------------------------------------------------------
 
-      use comcstfi_h, only: g, pi, r
+      use comcstfi_h, only: g, pi, r,rcp
       USE ioipsl_getin_p_mod, ONLY : getin_p
       use vertical_layers_mod, only : presnivs
@@ -62 +64 @@
     REAL, intent(in):: pv(ngrid,nlayer)  ! Meridional winds at full levels(m/s)
     REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) ! advected field nq
+    REAL, INTENT(IN) :: pht(ngrid,nlayer) ! advected field of potential temperature
+    REAL, INTENT(IN) :: pdht(ngrid,nlayer) ! tendancy of potential temperature
     REAL, INTENT(IN) :: pdq(ngrid,nlayer,nq)! tendancy field pq
     REAL,INTENT(in) :: pdt(ngrid,nlayer) ! Tendency on temperature (K/s)
@@ -74 +78 @@
     REAL, intent(out):: d_v(ngrid, nlayer) ! Tendency on meridional wind (m/s/s) due to gravity waves
     REAL, INTENT(out) :: d_pq(ngrid,nlayer,nq)! tendancy field pq
+    REAL :: d_h(ngrid, nlayer)  ! Tendency on PT (T/s/s) due to gravity waves mixing
     ! 0.3 INTERNAL ARRAYS
     REAL :: TT(ngrid, nlayer)   ! Temperature at full levels
@@ -79 +84 @@
     REAL :: UU(ngrid, nlayer)   ! Zonal wind at full levels
     REAL :: VV(ngrid, nlayer)   ! Meridional winds at full levels
+    REAL :: HH(ngrid, nlayer)   ! potential temperature at full levels
     REAL :: BVLOW(ngrid)        ! initial N at each grid (not used)
 
@@ -127 +133 @@
     REAL u_eddy_mix_p(NW, ngrid)       ! Zonal Diffusion coefficients
     REAL v_eddy_mix_p(NW, ngrid)       ! Meridional Diffusion coefficients
+    REAL h_eddy_mix_p(NW, ngrid)       ! potential temperature DC
     Real u_eddy_mix_tot(ngrid, nlayer+1)  ! Total zonal D
     Real v_eddy_mix_tot(ngrid, nlayer+1)  ! Total meridional D
+    Real h_eddy_mix_tot(ngrid, nlayer+1)  ! Total PT D
     REAL U_shear(ngrid,nlayer)
     Real wwp_vertical_tot(nlayer+1, NW, ngrid)  ! Total meridional D
@@ -246 +254 @@
     vv(:,:)=pv(:,:)+dtime*pdv(:,:)
     zq(:,:,:)=pq(:,:,:)+dtime*pdq(:,:,:)
+    hh(:,:)=pht(:,:)+dtime*pdht(:,:)
 ! Compute the real mass density by rho=p/R(T)T
      DO ll=1,nlayer
@@ -610 +619 @@
     u_eddy_mix_tot(:, :) = 0.
     v_eddy_mix_tot(:, :) = 0.
+    h_eddy_mix_tot(:, :) = 0.
     u_eddy_mix_p(:, :)=0.
     v_eddy_mix_p(:, :)=0.
+    h_eddy_mix_p(:, :)=0.
     d_eddy_mix_tot_ll(:,:,:)=0.
     pq_eddy_mix_p(:,:,:)=0.
@@ -627 +638 @@
                                /(ZH(:, LL + 1) - ZH(:, LL))                          &
                                *SIGN(1.,intr_freq_p(JW, :)) * SIN(ZP(JW, :))
+         h_eddy_mix_p(JW, :) = d_eddy_mix(JW,:)*(HH(:, LL + 1) - HH(:, LL))          &
+                               /(ZH(:, LL + 1) - ZH(:, LL))                          &
+                               *SIGN(1.,intr_freq_p(JW, :)) * SIN(ZP(JW, :))
+
       ENDDO
        u_eddy_mix_tot(:, LL+1)=0.
        v_eddy_mix_tot(:, LL+1)=0.
+       h_eddy_mix_tot(:, LL+1)=0.
       DO JW=1,NW
         u_eddy_mix_tot(:, LL+1) = u_eddy_mix_tot(:, LL+1) + u_eddy_mix_p(JW, :)
         v_eddy_mix_tot(:, LL+1) = v_eddy_mix_tot(:, LL+1) + v_eddy_mix_p(JW, :)
+        h_eddy_mix_tot(:, LL+1) = h_eddy_mix_tot(:, LL+1) + h_eddy_mix_p(JW, :)
       ENDDO
       DO JW=1,NW
@@ -699 +716 @@
     u_eddy_mix_tot(:, nlayer + 1) = 0.
     v_eddy_mix_tot(:, nlayer + 1) = 0.
+    h_eddy_mix_tot(:, nlayer + 1) = 0.
     pq_eddy_mix_tot(:, nlayer + 1,:)=0.
     ! Here, big change compared to FLott version:
@@ -706 +724 @@
       u_eddy_mix_tot(:, LL) = 0.
       v_eddy_mix_tot(:, LL) = 0.
+      h_eddy_mix_tot(:, LL) = 0.
     end DO
 
@@ -718 +737 @@
                               * (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
       v_eddy_mix_tot(:, LL) = v_eddy_mix_tot(:, LL - 1) + v_eddy_mix_tot(:, LAUNCH)     &
+                              * (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
+      h_eddy_mix_tot(:, LL) = h_eddy_mix_tot(:, LL - 1) + h_eddy_mix_tot(:, LAUNCH)     &
                               * (PH(:,LL)-PH(:,LL-1)) / (PH(:,LAUNCH)-PH(:,max(1,LAUNCH-3)))
     end DO
@@ -751 +772 @@
        !d_v(:, LL) =  (v_eddy_mix_tot(:, LL + 1) - v_eddy_mix_tot(:, LL)) &
        !             / (ZH(:, LL + 1) - ZH(:, LL)) 
+       d_h(:, LL) =  (h_eddy_mix_tot(:, LL + 1) - h_eddy_mix_tot(:, LL)) &
+                    / (ZH(:, LL + 1) - ZH(:, LL)) 
     ENDDO
 
@@ -761 +784 @@
     !df_eddymix_flx(:,:) = u_eddy_mix_tot(:,:)
     !d_pq(:, :, :)=0.
-    d_t(:,:) = 0.
+    !d_t(:,:) = 0.
     !d_v(:,:) = 0.
     !zustr(:) = 0.
@@ -779 +802 @@
     call write_output('du_eddymix_gwd','Tendency on U due to nonoro GW', 'm.s-2',du_eddymix_gwd(:,:))
     !call write_output('dv_eddymix_gwd','Tendency on V due to nonoro GW', 'm.s-2',dv_eddymix_gwd(:,:))   
-
-
+    d_h(:,:) = DTIME/DELTAT/REAL(NW) * d_h(:,:)                       &
+               + (1.-DTIME/DELTAT) * dh_eddymix_gwd(:,:)
+    do ii=1,ngrid
+    d_t(ii,:) = d_h(ii,:) * (PP(ii,:) / PH(ii,1))**rcp   
+    enddo                  
+               
+    dh_eddymix_gwd(:,:)=d_h(:,:)
+
+    DO QQ=1,NQ
+    d_pq(:, :, QQ) =DTIME/DELTAT/REAL(NW) * d_pq(:, :, QQ)            &
+               + (1.-DTIME/DELTAT) * dq_eddymix_gwd(:, :, QQ)
+    endDO
 
   END SUBROUTINE NONORO_GWD_MIX
@@ -789 +822 @@
 ! Subroutines used to allocate/deallocate module variables       
 ! ========================================================
-  SUBROUTINE ini_nonoro_gwd_mix(ngrid,nlayer)
+  SUBROUTINE ini_nonoro_gwd_mix(ngrid,nlayer,nq)
 
   IMPLICIT NONE
@@ -795 +828 @@
       INTEGER, INTENT (in) :: ngrid  ! number of atmospheric columns
       INTEGER, INTENT (in) :: nlayer ! number of atmospheric layers
+      INTEGER, INTENT (in) :: nq ! number of atmospheric tracers
 
          allocate(du_eddymix_gwd(ngrid,nlayer))
          allocate(dv_eddymix_gwd(ngrid,nlayer))
+         allocate(dh_eddymix_gwd(ngrid,nlayer))
+         allocate(dq_eddymix_gwd(ngrid,nlayer,nq))
          allocate(de_eddymix_rto(ngrid,nlayer+1))
          allocate(df_eddymix_flx(ngrid,nlayer+1))
@@ -816 +852 @@
     if (allocated(du_eddymix_gwd)) deallocate(du_eddymix_gwd)
     if (allocated(dv_eddymix_gwd)) deallocate(dv_eddymix_gwd)
+    if (allocated(dh_eddymix_gwd)) deallocate(dh_eddymix_gwd)
+    if (allocated(dq_eddymix_gwd)) deallocate(dq_eddymix_gwd)
     if (allocated(de_eddymix_rto)) deallocate(de_eddymix_rto)
     if (allocated(df_eddymix_flx)) deallocate(df_eddymix_flx)