Changeset 924

Timestamp:

06/19/19 17:22:32 (5 years ago)

Author:

dubos

Message:

devel : more options for the vertical profile of horizontal dissipation coefficient

File:

• codes/icosagcm/devel/src/dissip/dissip_gcm.f90 (modified) (5 diffs)

codes/icosagcm/devel/src/dissip/dissip_gcm.f90

@@ -26 +26 @@
   REAL,ALLOCATABLE,SAVE :: tau_divgrad(:)
 !$OMP THREADPRIVATE(tau_divgrad)
-  
+
   REAL,SAVE :: cgraddiv
 !$OMP THREADPRIVATE(cgraddiv)
@@ -33 +33 @@
   REAL,SAVE :: cdivgrad
 !$OMP THREADPRIVATE(cdivgrad)
-
+  
   INTEGER, SAVE :: rayleigh_friction_type, rayleigh_shear
 !$OMP THREADPRIVATE(rayleigh_friction_type)
@@ -126 +126 @@
        CALL dissip_constants
        CALL dissip_profile
+       CALL dissip_timescale
     END IF
 
@@ -236 +237 @@
   SUBROUTINE dissip_profile
     USE disvert_mod
-    INTEGER    :: l
-    REAL(rstd) :: mintau, zz, zvert, fact
-    fact=2
+    ! parameters used by the various profiles
+    ! IF planet_type == earth
+    !    IF dissip_vert_prof == 0
+    !      none
+    !    IF dissip_vert_prof == 1
+    !      dissip_zref, dissip_deltaz, dissip_factz
+    ! IF planet_type == other
+    !    IF dissip_vert_prof == 0
+    !      dissip_fac_mid
+    !      + dissip_deltaz, dissip_hdelta, dissip_fac_up, dissip_pupstart IF ok_strato     
+    !    IF dissip_vert_prof == 1
+    !      fac_mid, fac_up, startalt, delta => middle (hardcoded), scaleheight
+    !      
+
+    REAL(rstd), PARAMETER :: fact=2., &
+         fac_mid=3.,   & ! coefficient for lower/middle atmosphere
+         fac_up=30.,   & ! coefficient for upper atmosphere
+         startalt=70., & ! altitude (in km) for the transition from middle to upper atm.
+         delta=30.,    & ! Size (in km) of the transition region between middle/upper atm.
+         middle=startalt+delta/2.
+    REAL(rstd) :: dissip_zref, dissip_deltaz, dissip_factz, &
+         dissip_hdelta, dissip_fac_up, dissip_fac_mid, dissip_pupstart, &
+         scaleheight, &
+         zz, pseudoz, pup, &
+         sigma(llm), zvert(llm)
+    CHARACTER(LEN=255) :: planet_type ! earth, other, other_strato ; other_strato = other + ok_strato
+    LOGICAL :: ok_strato
+    INTEGER    :: l, dissip_vert_prof
+
+    ! select vertical profile of horizontal dissipation coefficients, see also [781]
+
+    planet_type='earth'
+    CALL getin('dissip_planet_type', planet_type)
+    SELECT CASE(planet_type)
+    CASE('earth','other')
+       ok_strato=.FALSE.
+    CASE('other_strato')
+       planet_type='other'
+       ok_strato=.TRUE.
+    CASE DEFAULT
+       STOP "Invalid value of dissip_planet_type, valid values are <earth>, <other>, <other_strato>"
+    END SELECT
+
+    dissip_vert_prof = 0
+    CALL getin('dissip_vert_prof',dissip_vert_prof)
+
+    SELECT CASE(dissip_vert_prof)
+    CASE(0)
+       IF(TRIM(planet_type)=='other') THEN
+          ! Default values given below are for a Venus-like planet
+          CALL getin('dissip_fac_mid',dissip_fac_mid )
+          dissip_fac_mid=2.
+          IF(ok_strato) THEN
+             dissip_fac_up=50.
+             ! deltaz et hdelta in km
+             dissip_deltaz=30.
+             dissip_hdelta=5.
+             ! pupstart in Pa
+             dissip_pupstart=1.e4          
+             CALL getin('dissip_deltaz',dissip_deltaz )
+             CALL getin('dissip_hdelta',dissip_hdelta )
+             CALL getin('dissip_fac_up',dissip_fac_up )
+             CALL getin('dissip_pupstart',dissip_pupstart)
+          END IF
+       END IF
+    CASE(1)
+       IF(TRIM(planet_type)=='earth') THEN
+          dissip_zref = 30.
+          dissip_deltaz = 10.
+          dissip_factz = 4.
+          CALL getin('dissip_zref',dissip_zref )
+          CALL getin('dissip_deltaz',dissip_deltaz )
+          CALL getin('dissip_factz',dissip_factz )
+       ELSE
+          ! fac_mid, fac_up, startalt, delta => middle are hardcoded
+          CALL getin('dissip_scaleheight', scaleheight)
+       END IF
+    CASE DEFAULT
+       STOP 'Invalid value of dissip_vert_prof : valid values are 0,1'
+    END SELECT
+
+    IF(ap_bp_present) THEN
+       sigma(:) = preff/presnivs(:)
+    ELSE
+       sigma(:) = 1.
+    END IF
+
+    SELECT CASE(TRIM(planet_type))
+    CASE('earth')
+       DO l=1,llm
+          IF(dissip_vert_prof == 1) THEN
+             pseudoz=8.*LOG(sigma(l))
+             zvert(l)=1+ &
+                  (TANH((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2. &
+                  *(dissip_factz-1.)
+          ELSE
+             zz = 1.-sigma(l)
+             zvert(l)= fact -( fact-1.)/( 1.+zz*zz )
+          END IF
+       END DO
+    CASE('other')
+       SELECT CASE(dissip_vert_prof)
+       CASE(0)
+          ! First step: going from 1 to dissip_fac_mid (in gcm.def)
+          !============
+          DO l=1,llm
+             zz       = 1. - sigma(l)
+             zvert(l) = dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
+          END DO
+          !
+          ! Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
+          !==========================
+          ! Utilisation de la fonction tangente hyperbolique pour augmenter
+          ! arbitrairement la dissipation et donc la stabilite du modele a 
+          ! partir d'une certaine altitude.
+          !
+          !   Le facteur multiplicatif de basse atmosphere etant deja pris 
+          !   en compte, il faut diviser le facteur multiplicatif de haute 
+          !   atmosphere par celui-ci.
+
+          IF(ok_strato) THEN
+             Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
+             DO l=1,llm
+                zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) &
+                     *(1.-(0.5*(1+tanh(-6./dissip_deltaz*              &
+                     (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
+             END DO
+          END IF
+       CASE(1)
+          DO l=1,llm
+             zz      = 1. - sigma(l)
+             zvert(l)= fac_mid -( fac_mid-1.)/( 1.+zz*zz )     
+             zvert(l)= zvert(l)*(1.0+((fac_up/fac_mid-1)*    &
+                  (1.-(0.5*(1+tanh(-6./                 &
+                  delta*(scaleheight*(-log(presnivs(l)/preff))-middle))))) &
+                  ))
+          END DO
+       END SELECT
+    END SELECT
+    
     DO l=1,llm
-       IF(ap_bp_present) THEN ! height-dependent dissipation
-          zz= 1. - preff/presnivs(l)
-       ELSE
-          zz = 0.
-       END IF
-       zvert=fact-(fact-1)/(1+zz*zz)
-       tau_graddiv(l) = tau_graddiv(l)/zvert
-       tau_gradrot(l) = tau_gradrot(l)/zvert
-       tau_divgrad(l) = tau_divgrad(l)/zvert
+       tau_graddiv(l) = tau_graddiv(l)/zvert(l)
+       tau_gradrot(l) = tau_gradrot(l)/zvert(l)
+       tau_divgrad(l) = tau_divgrad(l)/zvert(l)
     END DO
-    
+  END SUBROUTINE dissip_profile
+ 
+  SUBROUTINE dissip_timescale    
+    INTEGER :: l
+    REAL(rstd) :: mintau
     mintau=tau_graddiv(1)
     DO l=1,llm
@@ -276 +412 @@
     ENDIF
     
-  END SUBROUTINE dissip_profile
-  
+  END SUBROUTINE dissip_timescale
+
   SUBROUTINE dissip(f_ps,f_mass,f_phis,f_geopot,f_theta_rhodz,f_ue, f_dtheta_rhodz,f_due)
     TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_phis(:), f_geopot(:)