Changeset 1021 for trunk/LMDZ.COMMON


Ignore:
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
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F

    r1017 r1021  
    250250!Config  dissip_period>0 => on prend cette valeur
    251251       dissip_period = 0
     252       call getin('idissip',dissip_period) ! old Mars/Genreic model parameter
     253       ! if there is a "dissip_period" in run.def, it overrides "idissip"
    252254       CALL getin('dissip_period',dissip_period)
    253255
    … …  
    545547     $     "bad value for vert_prof_dissip")
    546548      else
    547        vert_prof_dissip=0
     549       vert_prof_dissip=0 ! default for planets !
     550       if (planet_type=="mars") then
     551         vert_prof_dissip=1 ! use fac_mid & fac_up & startalt & delta
     552       endif
    548553      endif
    549554

  • trunk/LMDZ.COMMON/libf/dyn3d/inidissip.F90

    r979 r1021  
    55     tetagdiv,tetagrot,tetatemp, vert_prof_dissip)
    66  !=======================================================================
    7   !   initialisation de la dissipation horizontale
     7  !   Initialization for horizontal (lateral) dissipation
     8  !  - in all cases, there is a multiplicative coefficient which increases
     9  !    the dissipation in the higher levels of the atmosphere, but there
     10  !    are different ways of seting the vertical profile of this coefficient
     11  !    (see code below).
     12  !  - the call to dissipation, every 'dissip_period' dynamical time step,
     13  !    can be imposed via 'dissip_period=...' in run.def (or via the
     14  !    'idissip=...' flag, but this flag should become obsolete, and is
     15  !    overridden by the 'dissip_period' flag). Note that setting
     16  !    dissip_period=0 has the special meaning of requesting an "optimal"
     17  !    value for "dissip_period" (then taken as the largest possible value)
     18  !  - the three characteristic time scales (relative to the smallest
     19  !    longitudinal grid mesh size), which are privided in run.def, which
     20  !    are used for the dissipations steps are:
     21  !     tetagdiv : time scale for the gradient of the divergence of velocities
     22  !     tetagrot : time scale for the curl of the curl of velocities
     23  !     tetatemp : time scale for the laplacian of the potential temperature
    824  !=======================================================================
    9   !-----------------------------------------------------------------------
    10   !   declarations:
    11   !   -------------
    1225
    1326  USE control_mod, only : dissip_period,iperiod,planet_type
    … …  
    2639  REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp
    2740
    28   integer, INTENT(in):: vert_prof_dissip ! for the Earth model !!
    29   ! Vertical profile of horizontal dissipation
     41  integer, INTENT(in):: vert_prof_dissip ! Vertical profile of horizontal dissipation
     42  ! For the Earth model:
    3043  ! Allowed values:
    3144  ! 0: rational fraction, function of pressure
    3245  ! 1: tanh of altitude
    33 
     46  ! For planets:
     47  ! 0: use fac_mid (read from run.def)
     48  ! 1: use fac_mid, fac_up, startalt, delta (hard coded in inidissip.F90)
    3449! Local variables:
    3550  REAL fact,zvert(llm),zz
    … …  
    4661
    4762  REAL ran1
    48 
     63  logical,save :: firstcall=.true.
     64  real,save :: fac_mid,fac_up,startalt,delta,middle
     65
     66  if (firstcall) then
     67    firstcall=.false.
     68    if ((planet_type.ne."earth").and.(vert_prof_dissip.eq.1)) then
     69      ! initialize values for dissipation variation along the vertical (Mars)
     70      fac_mid=3 ! coefficient for lower/middle atmosphere
     71      fac_up=30 ! coefficient for upper atmosphere
     72      startalt=70. ! altitude (in km) for the transition from middle to upper atm.
     73      delta=30.! Size (in km) of the transition region between middle/upper atm.
     74
     75      middle=startalt+delta/2
     76      write(lunout,*)"inidissip: multiplicative factors in altitude:", &
     77        " fac_mid=",fac_mid," fac_up=",fac_up
     78      write(lunout,*)" transition mid/up : startalt (km) =",startalt, &
     79        " delta (km) =",delta
     80    endif
     81  endif !of if firstcall
    4982
    5083  !-----------------------------------------------------------------------
    … …  
    196229  else ! other planets
    197230 
     231   if (vert_prof_dissip == 0) then
    198232! First step: going from 1 to dissip_fac_mid (in gcm.def)
    199233!============
    200    DO l=1,llm
     234    DO l=1,llm
    201235     zz      = 1. - preff/presnivs(l)
    202236     zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
    203    ENDDO
    204 
    205    write(lunout,*) 'Dissipation : '
    206    write(lunout,*) 'Multiplication de la dissipation en altitude :'
    207    write(lunout,*) '  dissip_fac_mid =', dissip_fac_mid
     237    ENDDO
     238
     239    write(lunout,*) 'Dissipation : '
     240    write(lunout,*) 'Multiplication de la dissipation en altitude :'
     241    write(lunout,*) '  dissip_fac_mid =', dissip_fac_mid
    208242
    209243! Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
    … …  
    217251!   atmosphere par celui-ci.
    218252
    219    if (ok_strato) then
    220 
    221     Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
    222     do l=1,llm
     253    if (ok_strato) then
     254
     255     Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
     256     do l=1,llm
    223257      zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) &
    224258                *(1-(0.5*(1+tanh(-6./dissip_deltaz*              &
    225259               (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
    226     enddo
    227 
    228     write(*,*) '  dissip_fac_up =', dissip_fac_up
    229     write(*,*) 'Transition mid /up:  Pupstart,delta =',           &
     260     enddo
     261
     262     write(*,*) '  dissip_fac_up =', dissip_fac_up
     263     write(*,*) 'Transition mid /up:  Pupstart,delta =',           &
    230264                   dissip_pupstart,'Pa', dissip_deltaz , '(km)'
    231265
    232    endif ! of if (ok_strato)
    233  
     266    endif ! of if (ok_strato)
     267   elseif (vert_prof_dissip==1) then
     268    DO l=1,llm
     269     zz      = 1. - preff/presnivs(l)
     270     zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
     271     
     272     zvert(l)= zvert(l)*(1.0+((fac_up/fac_mid-1)*    &
     273                (1-(0.5*(1+tanh(-6./                 &
     274                delta*(scaleheight*(-log(presnivs(l)/preff))-middle))))) &
     275                ))
     276    ENDDO
     277
     278   else
     279     write(lunout,*) 'wrong value for vert_prof_dissip:',vert_prof_dissip
     280     abort_message='wrong value for vert_prof_dissip'
     281     call abort_gcm('inidissip',abort_message,1)
     282   endif ! of (vert_prof_dissip == 0)
    234283  endif ! of if (planet_type.eq."earth")
    235284
    236285
    237   write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale'
     286  write(lunout,*)'inidissip: Time constants for lateral dissipation'
    238287
    239288  tetamin =  1.e+6
    … …  
    259308  write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr
    260309
     310  write(lunout,*)'pseudoZ(km)  zvert    dt(tetagdiv)   dt(tetagrot)   dt(divgrad)'
    261311  DO l = 1,llm
    262      write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), &
    263           dtdiss*tetah(l)
    264   ENDDO
     312     write(lunout,'(f6.1,x,4(1pe14.7,x))') &
     313     pseudoalt(l),zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l),dtdiss*tetah(l)
     314     ! test if disipation is not too strong (for Explicit Euler time marching)
     315     if (dtdiss*tetaudiv(l).gt.1.9) then
     316       write(lunout,*)"STOP : lateral dissipation is too intense and will"
     317       write(lunout,*)"       generate instabilities in the model !"
     318       write(lunout,*)" You must increase tetagdiv (or increase dissip_period"
     319       write(lunout,*)"                             or increase day_stap)"
     320     endif
     321     if (dtdiss*tetaurot(l).gt.1.9) then
     322       write(lunout,*)"STOP : lateral dissipation is too intense and will"
     323       write(lunout,*)"       generate instabilities in the model !"
     324       write(lunout,*)" You must increase tetaurot (or increase dissip_period"
     325       write(lunout,*)"                             or increase day_stap)"
     326     endif
     327     if (dtdiss*tetah(l).gt.1.9) then
     328       write(lunout,*)"STOP : lateral dissipation is too intense and will"
     329       write(lunout,*)"       generate instabilities in the model !"
     330       write(lunout,*)" You must increase tetah (or increase dissip_period"
     331       write(lunout,*)"                          or increase day_stap)"
     332     endif
     333  ENDDO ! of DO l=1,llm
    265334
    266335END SUBROUTINE inidissip

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

    r1019 r1021  
    277277!Config  dissip_period>0 => on prend cette valeur
    278278       dissip_period = 0
     279       call getin('idissip',dissip_period) ! old Mars/Genreic model parameter
     280       ! if there is a "dissip_period" in run.def, it overrides "idissip"
    279281       CALL getin('dissip_period',dissip_period)
    280282
    … …  
    601603     $     "bad value for vert_prof_dissip")
    602604      else
    603        vert_prof_dissip=0
     605       vert_prof_dissip=0 ! default for planets !
     606       if (planet_type=="mars") then
     607         vert_prof_dissip=1 ! use fac_mid & fac_up & startalt & delta
     608       endif
    604609      endif
    605610

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

    r979 r1021  
    55     tetagdiv,tetagrot,tetatemp, vert_prof_dissip)
    66  !=======================================================================
    7   !   initialisation de la dissipation horizontale
     7  !   Initialization for horizontal (lateral) dissipation
     8  !  - in all cases, there is a multiplicative coefficient which increases
     9  !    the dissipation in the higher levels of the atmosphere, but there
     10  !    are different ways of seting the vertical profile of this coefficient
     11  !    (see code below).
     12  !  - the call to dissipation, every 'dissip_period' dynamical time step,
     13  !    can be imposed via 'dissip_period=...' in run.def (or via the
     14  !    'idissip=...' flag, but this flag should become obsolete, and is
     15  !    overridden by the 'dissip_period' flag). Note that setting
     16  !    dissip_period=0 has the special meaning of requesting an "optimal"
     17  !    value for "dissip_period" (then taken as the largest possible value)
     18  !  - the three characteristic time scales (relative to the smallest
     19  !    longitudinal grid mesh size), which are privided in run.def, which
     20  !    are used for the dissipations steps are:
     21  !     tetagdiv : time scale for the gradient of the divergence of velocities
     22  !     tetagrot : time scale for the curl of the curl of velocities
     23  !     tetatemp : time scale for the laplacian of the potential temperature
    824  !=======================================================================
    9   !-----------------------------------------------------------------------
    10   !   declarations:
    11   !   -------------
    1225
    1326  USE control_mod, only : dissip_period,iperiod,planet_type
    … …  
    2639  REAL,INTENT(in) :: tetagdiv,tetagrot,tetatemp
    2740
    28   integer, INTENT(in):: vert_prof_dissip ! for the Earth model !!
    29   ! Vertical profile of horizontal dissipation
     41  integer, INTENT(in):: vert_prof_dissip ! Vertical profile of horizontal dissipation
     42  ! For the Earth model:
    3043  ! Allowed values:
    3144  ! 0: rational fraction, function of pressure
    3245  ! 1: tanh of altitude
    33 
     46  ! For planets:
     47  ! 0: use fac_mid (read from run.def)
     48  ! 1: use fac_mid, fac_up, startalt, delta (hard coded in inidissip.F90)
    3449! Local variables:
    3550  REAL fact,zvert(llm),zz
    … …  
    4661
    4762  REAL ran1
    48 
     63  logical,save :: firstcall=.true.
     64  real,save :: fac_mid,fac_up,startalt,delta,middle
     65
     66  if (firstcall) then
     67    firstcall=.false.
     68    if ((planet_type.ne."earth").and.(vert_prof_dissip.eq.1)) then
     69      ! initialize values for dissipation variation along the vertical (Mars)
     70      fac_mid=3 ! coefficient for lower/middle atmosphere
     71      fac_up=30 ! coefficient for upper atmosphere
     72      startalt=70. ! altitude (in km) for the transition from middle to upper atm.
     73      delta=30.! Size (in km) of the transition region between middle/upper atm.
     74
     75      middle=startalt+delta/2
     76      write(lunout,*)"inidissip: multiplicative factors in altitude:", &
     77        " fac_mid=",fac_mid," fac_up=",fac_up
     78      write(lunout,*)" transition mid/up : startalt (km) =",startalt, &
     79        " delta (km) =",delta
     80    endif
     81  endif !of if firstcall
    4982
    5083  !-----------------------------------------------------------------------
    … …  
    196229  else ! other planets
    197230 
     231   if (vert_prof_dissip == 0) then
    198232! First step: going from 1 to dissip_fac_mid (in gcm.def)
    199233!============
    200    DO l=1,llm
     234    DO l=1,llm
    201235     zz      = 1. - preff/presnivs(l)
    202236     zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
    203    ENDDO
    204 
    205    write(lunout,*) 'Dissipation : '
    206    write(lunout,*) 'Multiplication de la dissipation en altitude :'
    207    write(lunout,*) '  dissip_fac_mid =', dissip_fac_mid
     237    ENDDO
     238
     239    write(lunout,*) 'Dissipation : '
     240    write(lunout,*) 'Multiplication de la dissipation en altitude :'
     241    write(lunout,*) '  dissip_fac_mid =', dissip_fac_mid
    208242
    209243! Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
    … …  
    217251!   atmosphere par celui-ci.
    218252
    219    if (ok_strato) then
    220 
    221     Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
    222     do l=1,llm
     253    if (ok_strato) then
     254
     255     Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
     256     do l=1,llm
    223257      zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1) &
    224258                *(1-(0.5*(1+tanh(-6./dissip_deltaz*              &
    225259               (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
    226     enddo
    227 
    228     write(*,*) '  dissip_fac_up =', dissip_fac_up
    229     write(*,*) 'Transition mid /up:  Pupstart,delta =',           &
     260     enddo
     261
     262     write(*,*) '  dissip_fac_up =', dissip_fac_up
     263     write(*,*) 'Transition mid /up:  Pupstart,delta =',           &
    230264                   dissip_pupstart,'Pa', dissip_deltaz , '(km)'
    231265
    232    endif ! of if (ok_strato)
    233  
     266    endif ! of if (ok_strato)
     267   elseif (vert_prof_dissip==1) then
     268    DO l=1,llm
     269     zz      = 1. - preff/presnivs(l)
     270     zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
     271     
     272     zvert(l)= zvert(l)*(1.0+((fac_up/fac_mid-1)*    &
     273                (1-(0.5*(1+tanh(-6./                 &
     274                delta*(scaleheight*(-log(presnivs(l)/preff))-middle))))) &
     275                ))
     276    ENDDO
     277
     278   else
     279     write(lunout,*) 'wrong value for vert_prof_dissip:',vert_prof_dissip
     280     abort_message='wrong value for vert_prof_dissip'
     281     call abort_gcm('inidissip',abort_message,1)
     282   endif ! of (vert_prof_dissip == 0)
    234283  endif ! of if (planet_type.eq."earth")
    235284
    236285
    237   write(lunout,*)'inidissip: Constantes de temps de la diffusion horizontale'
     286  write(lunout,*)'inidissip: Time constants for lateral dissipation'
    238287
    239288  tetamin =  1.e+6
    … …  
    259308  write(lunout,*)'inidissip: dissip_period=',dissip_period,' dtdiss=',dtdiss,' dtvr=',dtvr
    260309
     310  write(lunout,*)'pseudoZ(km)  zvert    dt(tetagdiv)   dt(tetagrot)   dt(divgrad)'
    261311  DO l = 1,llm
    262      write(lunout,*)zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l), &
    263           dtdiss*tetah(l)
    264   ENDDO
     312     write(lunout,'(f6.1,x,4(1pe14.7,x))') &
     313     pseudoalt(l),zvert(l),dtdiss*tetaudiv(l),dtdiss*tetaurot(l),dtdiss*tetah(l)
     314     ! test if disipation is not too strong (for Explicit Euler time marching)
     315     if (dtdiss*tetaudiv(l).gt.1.9) then
     316       write(lunout,*)"STOP : lateral dissipation is too intense and will"
     317       write(lunout,*)"       generate instabilities in the model !"
     318       write(lunout,*)" You must increase tetagdiv (or increase dissip_period"
     319       write(lunout,*)"                             or increase day_stap)"
     320     endif
     321     if (dtdiss*tetaurot(l).gt.1.9) then
     322       write(lunout,*)"STOP : lateral dissipation is too intense and will"
     323       write(lunout,*)"       generate instabilities in the model !"
     324       write(lunout,*)" You must increase tetaurot (or increase dissip_period"
     325       write(lunout,*)"                             or increase day_stap)"
     326     endif
     327     if (dtdiss*tetah(l).gt.1.9) then
     328       write(lunout,*)"STOP : lateral dissipation is too intense and will"
     329       write(lunout,*)"       generate instabilities in the model !"
     330       write(lunout,*)" You must increase tetah (or increase dissip_period"
     331       write(lunout,*)"                          or increase day_stap)"
     332     endif
     333  ENDDO ! of DO l=1,llm
    265334
    266335END SUBROUTINE inidissip
Note: See TracChangeset for help on using the changeset viewer.