Changeset 1021 for trunk/LMDZ.COMMON/libf/dyn3dpar

Timestamp:

Aug 28, 2013, 9:19:53 AM (11 years ago)

Author:

emillour

Message:

Common dynamics: add the possibility to handle dissipation as in the Mars GCM.
EM

Location:

trunk/LMDZ.COMMON/libf/dyn3dpar

Files:

• conf_gcm.F (modified) (2 diffs)
• inidissip.F90 (modified) (6 diffs)

trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F

@@ -277 +277 @@
 !Config  dissip_period>0 => on prend cette valeur
        dissip_period = 0
+       call getin('idissip',dissip_period) ! old Mars/Genreic model parameter
+       ! if there is a "dissip_period" in run.def, it overrides "idissip"
        CALL getin('dissip_period',dissip_period)
 
@@ -601 +603 @@
      $     "bad value for vert_prof_dissip")
       else
-       vert_prof_dissip=0
+       vert_prof_dissip=0 ! default for planets !
+       if (planet_type=="mars") then
+         vert_prof_dissip=1 ! use fac_mid & fac_up & startalt & delta
+       endif
       endif
 

trunk/LMDZ.COMMON/libf/dyn3dpar/inidissip.F90

@@ -5 +5 @@
      tetagdiv,tetagrot,tetatemp, vert_prof_dissip)
   !=======================================================================
-  !   initialisation de la dissipation horizontale
+  !   Initialization for horizontal (lateral) dissipation
+  !  - in all cases, there is a multiplicative coefficient which increases
+  !    the dissipation in the higher levels of the atmosphere, but there
+  !    are different ways of seting the vertical profile of this coefficient
+  !    (see code below).
+  !  - the call to dissipation, every 'dissip_period' dynamical time step,
+  !    can be imposed via 'dissip_period=...' in run.def (or via the
+  !    'idissip=...' flag, but this flag should become obsolete, and is
+  !    overridden by the 'dissip_period' flag). Note that setting 
+  !    dissip_period=0 has the special meaning of requesting an "optimal"
+  !    value for "dissip_period" (then taken as the largest possible value)
+  !  - the three characteristic time scales (relative to the smallest
+  !    longitudinal grid mesh size), which are privided in run.def, which 
+  !    are used for the dissipations steps are:
+  !     tetagdiv : time scale for the gradient of the divergence of velocities
+  !     tetagrot : time scale for the curl of the curl of velocities
+  !     tetatemp : time scale for the laplacian of the potential temperature
   !=======================================================================
-  !-----------------------------------------------------------------------
-  !   declarations:
-  !   -------------
 
   USE control_mod, only : dissip_period,iperiod,planet_type
@@ -26 +39 @@
   REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp
 
-  integer, INTENT(in):: vert_prof_dissip ! for the Earth model !!
-  ! Vertical profile of horizontal dissipation
+  integer, INTENT(in):: vert_prof_dissip ! Vertical profile of horizontal dissipation
+  ! For the Earth model:
   ! Allowed values:
   ! 0: rational fraction, function of pressure
   ! 1: tanh of altitude
-
+  ! For planets:
+  ! 0: use fac_mid (read from run.def)
+  ! 1: use fac_mid, fac_up, startalt, delta (hard coded in inidissip.F90)
 ! Local variables:
   REAL fact,zvert(llm),zz
@@ -46 +61 @@
 
   REAL ran1
-
+  logical,save :: firstcall=.true.
+  real,save :: fac_mid,fac_up,startalt,delta,middle
+
+  if (firstcall) then
+    firstcall=.false.
+    if ((planet_type.ne."earth").and.(vert_prof_dissip.eq.1)) then
+      ! initialize values for dissipation variation along the vertical (Mars)
+      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
+      write(lunout,*)"inidissip: multiplicative factors in altitude:", &
+        " fac_mid=",fac_mid," fac_up=",fac_up
+      write(lunout,*)" transition mid/up : startalt (km) =",startalt, &
+        " delta (km) =",delta
+    endif
+  endif !of if firstcall
 
   !-----------------------------------------------------------------------
@@ -196 +229 @@
   else ! other planets
   
+   if (vert_prof_dissip == 0) then
 ! First step: going from 1 to dissip_fac_mid (in gcm.def)
 !============
-   DO l=1,llm
+    DO l=1,llm
      zz      = 1. - preff/presnivs(l)
      zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
-   ENDDO
-
-   write(lunout,*) 'Dissipation : '
-   write(lunout,*) 'Multiplication de la dissipation en altitude :'
-   write(lunout,*) '  dissip_fac_mid =', dissip_fac_mid
+    ENDDO
+
+    write(lunout,*) 'Dissipation : '
+    write(lunout,*) 'Multiplication de la dissipation en altitude :'
+    write(lunout,*) '  dissip_fac_mid =', dissip_fac_mid
 
 ! Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
@@ -217 +251 @@
 !   atmosphere par celui-ci.
 
-   if (ok_strato) then
-
-    Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
-    do l=1,llm
+    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))  ))))  ))
-    enddo 
-
-    write(*,*) '  dissip_fac_up =', dissip_fac_up
-    write(*,*) 'Transition mid /up:  Pupstart,delta =',           &
+     enddo 
+
+     write(*,*) '  dissip_fac_up =', dissip_fac_up
+     write(*,*) 'Transition mid /up:  Pupstart,delta =',           &
                    dissip_pupstart,'Pa', dissip_deltaz , '(km)'
 
-   endif ! of if (ok_strato)
-  
+    endif ! of if (ok_strato)
+   elseif (vert_prof_dissip==1) then
+    DO l=1,llm
+     zz      = 1. - preff/presnivs(l)
+     zvert(l)= dissip_fac_mid -( dissip_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))))) &
+                ))
+    ENDDO
+
+   else
+     write(lunout,*) 'wrong value for vert_prof_dissip:',vert_prof_dissip
+     abort_message='wrong value for vert_prof_dissip'
+     call abort_gcm('inidissip',abort_message,1)
+   endif ! of (vert_prof_dissip == 0)
   endif ! of if (planet_type.eq."earth")
 
 
-  write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale'
+  write(lunout,*)'inidissip: Time constants for lateral dissipation'
 
   tetamin =  1.e+6
@@ -259 +308 @@
   write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr
 
+  write(lunout,*)'pseudoZ(km)  zvert    dt(tetagdiv)   dt(tetagrot)   dt(divgrad)'
   DO l = 1,llm
-     write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), &
-          dtdiss*tetah(l)
-  ENDDO
+     write(lunout,'(f6.1,x,4(1pe14.7,x))') &
+     pseudoalt(l),zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l),dtdiss*tetah(l)
+     ! test if disipation is not too strong (for Explicit Euler time marching)
+     if (dtdiss*tetaudiv(l).gt.1.9) then
+       write(lunout,*)"STOP : lateral dissipation is too intense and will"
+       write(lunout,*)"       generate instabilities in the model !"
+       write(lunout,*)" You must increase tetagdiv (or increase dissip_period"
+       write(lunout,*)"                             or increase day_stap)"
+     endif
+     if (dtdiss*tetaurot(l).gt.1.9) then
+       write(lunout,*)"STOP : lateral dissipation is too intense and will"
+       write(lunout,*)"       generate instabilities in the model !"
+       write(lunout,*)" You must increase tetaurot (or increase dissip_period"
+       write(lunout,*)"                             or increase day_stap)"
+     endif
+     if (dtdiss*tetah(l).gt.1.9) then
+       write(lunout,*)"STOP : lateral dissipation is too intense and will"
+       write(lunout,*)"       generate instabilities in the model !"
+       write(lunout,*)" You must increase tetah (or increase dissip_period"
+       write(lunout,*)"                          or increase day_stap)"
+     endif
+  ENDDO ! of DO l=1,llm
 
 END SUBROUTINE inidissip