Ignore:
Timestamp:
Sep 30, 2010, 10:29:10 AM (14 years ago)
Author:
Ehouarn Millour
Message:

Implementation of FC's improved Newtonian configuration: more realistic (and tunable) friction & fields towards which to relax and planet parameters (gravity, radius, rotation rate...) can be read in from a planet.def file (this occurs only when running in "newtonian mode"; planet parameters are otherwise read in from the start file).

Checked that these implementation don't change results when running standard Earth bench tests on Vargas.

EM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/iniacademic.F

    r1403 r1437  
    99      USE infotrac, ONLY : nqtot
    1010      USE control_mod
     11#ifdef CPP_IOIPSL
     12      USE IOIPSL
     13#else
     14! if not using IOIPSL, we still need to use (a local version of) getin
     15      USE ioipsl_getincom
     16#endif
     17      USE Write_Field
    1118 
    1219
    … …  
    7077      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
    7178      REAL phi(ip1jmp1,llm)                  ! geopotentiel
    72       REAL ddsin,tetarappelj,tetarappell,zsig
     79      REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
    7380      real tetajl(jjp1,llm)
    7481      INTEGER i,j,l,lsup,ij
    7582
     83      REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
     84      REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
     85      LOGICAL ok_geost             ! Initialisation vent geost. ou nul
     86      LOGICAL ok_pv                ! Polar Vortex
     87      REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex
     88     
    7689      real zz,ran1
    7790      integer idum
    … …  
    8497      if (planet_type=="earth") then
    8598c
     99        ! initialize planet radius, rotation rate,...
     100        call conf_planete
     101
    86102        time_0=0.
    87         day_ref=0
     103        day_ref=1
    88104        annee_ref=0
    89105
    90106        im         = iim
    91107        jm         = jjm
    92         day_ini    = 0
    93         omeg       = 4.*asin(1.)/86400.
    94         rad    = 6371229.
    95         g      = 9.8
    96         daysec = 86400.
     108        day_ini    = 1
    97109        dtvr    = daysec/REAL(day_step)
    98110        zdtvr=dtvr
    99         kappa  = 0.2857143
    100         cpp    = 1004.70885
    101         preff     = 101325.
    102         pa        =  50000.
    103111        etot0      = 0.
    104112        ptot0      = 0.
    … …  
    129137        if (iflag_phys.eq.2) then
    130138          ! initializations for the academic case
    131           ps(:)=1.e5
    132           phis(:)=0.
    133 c---------------------------------------------------------------------
    134 
    135           taurappel=10.*daysec
    136 
    137 c---------------------------------------------------------------------
    138 c   Calcul de la temperature potentielle :
    139 c   --------------------------------------
    140 
     139         
     140          ! 1. local parameters
     141          ! by convention, winter is in the southern hemisphere
     142          ! Geostrophic wind or no wind?
     143          ok_geost=.TRUE.
     144          CALL getin('ok_geost',ok_geost)
     145          ! Constants for Newtonian relaxation and friction
     146          k_f=1.                !friction
     147          CALL getin('k_j',k_f)
     148          k_f=1./(daysec*k_f)
     149          k_c_s=4.  !cooling surface
     150          CALL getin('k_c_s',k_c_s)
     151          k_c_s=1./(daysec*k_c_s)
     152          k_c_a=40. !cooling free atm
     153          CALL getin('k_c_a',k_c_a)
     154          k_c_a=1./(daysec*k_c_a)
     155          ! Constants for Teta equilibrium profile
     156          teta0=315.     ! mean Teta (S.H. 315K)
     157          CALL getin('teta0',teta0)
     158          ttp=200.       ! Tropopause temperature (S.H. 200K)
     159          CALL getin('ttp',ttp)
     160          eps=0.         ! Deviation to N-S symmetry(~0-20K)
     161          CALL getin('eps',eps)
     162          delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
     163          CALL getin('delt_y',delt_y)
     164          delt_z=10.     ! Vertical Gradient (S.H. 10K)
     165          CALL getin('delt_z',delt_z)
     166          ! Polar vortex
     167          ok_pv=.false.
     168          CALL getin('ok_pv',ok_pv)
     169          phi_pv=-50.            ! Latitude of edge of vortex
     170          CALL getin('phi_pv',phi_pv)
     171          phi_pv=phi_pv*pi/180.
     172          dphi_pv=5.             ! Width of the edge
     173          CALL getin('dphi_pv',dphi_pv)
     174          dphi_pv=dphi_pv*pi/180.
     175          gam_pv=4.              ! -dT/dz vortex (in K/km)
     176          CALL getin('gam_pv',gam_pv)
     177         
     178          ! 2. Initialize fields towards which to relax
     179          ! Friction
     180          knewt_g=k_c_a
    141181          DO l=1,llm
    142             zsig=ap(l)/preff+bp(l)
    143             if (zsig.gt.0.3) then
    144              lsup=l
    145              tetarappell=1./8.*(-log(zsig)-.5)
    146              DO j=1,jjp1
    147              ddsin=sin(rlatu(j))-sin(pi/20.)
    148              tetajl(j,l)=300.*(1+1./18.*(1.-3.*ddsin*ddsin)+tetarappell)
    149              ENDDO
    150             else
    151 c   Choix isotherme au-dessus de 300 mbar
    152              do j=1,jjp1
    153                tetajl(j,l)=tetajl(j,lsup)*(0.3/zsig)**kappa
    154              enddo
    155             endif ! of if (zsig.gt.0.3)
     182           zsig=presnivs(l)/preff
     183           knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
     184           kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
     185          ENDDO
     186          DO j=1,jjp1
     187            clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
     188          ENDDO
     189         
     190          ! Potential temperature
     191          DO l=1,llm
     192           zsig=presnivs(l)/preff
     193           tetastrat=ttp*zsig**(-kappa)
     194           tetapv=tetastrat
     195           IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
     196               tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
     197           ENDIF
     198           DO j=1,jjp1
     199             ! Troposphere
     200             ddsin=sin(rlatu(j))
     201             tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin
     202     &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
     203             ! Profil stratospherique isotherme (+vortex)
     204             w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
     205             tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
     206             tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
     207           ENDDO
    156208          ENDDO ! of DO l=1,llm
     209 
     210!          CALL writefield('theta_eq',tetajl)
    157211
    158212          do l=1,llm
    … …  
    165219          enddo
    166220
    167 c       call dump2d(jjp1,llm,tetajl,'TEQ   ')
    168 
    169           CALL pression ( ip1jmp1, ap, bp, ps, p       )
    170           CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    171           CALL massdair(p,masse)
    172 
    173 c  intialisation du vent et de la temperature
    174           teta(:,:)=tetarappel(:,:)
    175           CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
    176           call ugeostr(phi,ucov)
    177           vcov=0.
    178           q(:,:,1   )=1.e-10
    179           q(:,:,2   )=1.e-15
    180           q(:,:,3:nqtot)=0.
    181 
    182 
    183 c   perturbation aleatoire sur la temperature
    184           idum  = -1
    185           zz = ran1(idum)
    186           idum  = 0
    187           do l=1,llm
    188             do ij=iip2,ip1jm
    189               teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
    190             enddo
    191           enddo
    192 
    193           do l=1,llm
    194             do ij=1,ip1jmp1,iip1
    195               teta(ij+iim,l)=teta(ij,l)
    196             enddo
    197           enddo
    198 
    199 
     221          ! 3. Initialize fields (if necessary)
     222          IF (.NOT. read_start) THEN
     223            ! surface pressure
     224            ps(:)=preff
     225            ! ground geopotential
     226            phis(:)=0.
     227           
     228            CALL pression ( ip1jmp1, ap, bp, ps, p       )
     229            CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     230            CALL massdair(p,masse)
     231
     232            ! bulk initialization of temperature
     233            teta(:,:)=tetarappel(:,:)
     234           
     235            ! geopotential
     236            CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
     237           
     238            ! winds
     239            if (ok_geost) then
     240              call ugeostr(phi,ucov)
     241            else
     242              ucov(:,:)=0.
     243            endif
     244            vcov(:,:)=0.
     245           
     246            ! bulk initialization of tracers
     247            do i=1,nqtot
     248              if (i.eq.1) q(:,:,i)=1.e-10
     249              if (i.eq.2) q(:,:,i)=1.e-15
     250              if (i.gt.2) q(:,:,i)=0.
     251            enddo
     252
     253            ! add random perturbation to temperature
     254            idum  = -1
     255            zz = ran1(idum)
     256            idum  = 0
     257            do l=1,llm
     258              do ij=iip2,ip1jm
     259                teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
     260              enddo
     261            enddo
     262
     263            do l=1,llm
     264              do ij=1,ip1jmp1,iip1
     265                teta(ij+iim,l)=teta(ij,l)
     266              enddo
     267            enddo
    200268
    201269c     PRINT *,' Appel test_period avec tetarappel '
    … …  
    204272c     CALL  test_period ( ucov,vcov,teta,q,p,phis )
    205273
    206 c   initialisation d'un traceur sur une colonne
    207           j=jjp1*3/4
    208           i=iip1/2
    209           ij=(j-1)*iip1+i
    210           q(ij,:,3)=1.
     274           ! initialize a traceur on one column
     275!          j=jjp1*3/4
     276!          i=iip1/2
     277!          ij=(j-1)*iip1+i
     278!          q(ij,:,3)=1.
     279
     280          ENDIF ! of IF (.NOT. read_start)
    211281        endif ! of if (iflag_phys.eq.2)
    212282       
Note: See TracChangeset for help on using the changeset viewer.