Changeset 1437 for LMDZ5/branches


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

Location:
LMDZ5/branches/LMDZ5V1.0-dev/libf
Files:
2 added
10 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/academic.h

    r524 r1437  
    11!
    2 ! $Header$
     2! $Id$
    33!
    4       real tetarappel(ip1jmp1,llm),taurappel
    5       common/academic/tetarappel,taurappel
     4      common/academic/tetarappel,knewt_t,kfrict,knewt_g,clat4
     5      real :: tetarappel(ip1jmp1,llm)
     6      real :: knewt_t(llm)
     7      real :: kfrict(llm)
     8      real :: knewt_g
     9      real :: clat4(ip1jmp1)

  • LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/comconst.h

    r1279 r1437  
    55! INCLUDE comconst.h
    66
    7       COMMON/comconst/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,           &
    8      & dtvr,daysec,                                                     &
     7      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
     8     &                 iflag_top_bound
     9      COMMON/comconstr/dtvr,daysec,                                     &
    910     & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg              &
    1011     &                   ,dissip_factz,dissip_deltaz,dissip_zref        &
    11      &                   ,iflag_top_bound,tau_top_bound
     12     &                   ,tau_top_bound,                                &
     13     & daylen,year_day,molmass
    1214
    1315
    1416      INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
    15       REAL dtvr,daysec
    16       REAL pi,dtphys,dtdiss,rad,r,cpp,kappa
    17       REAL cotot,unsim,g,omeg
     17      REAL dtvr ! dynamical time step (in s)
     18      REAL daysec !length (in s) of a standard day
     19      REAL pi    ! something like 3.14159....
     20      REAL dtphys ! (s) time step for the physics
     21      REAL dtdiss ! (s) time step for the dissipation
     22      REAL rad ! (m) radius of the planet
     23      REAL r ! Gas constant R=8.31 J.K-1.mol-1
     24      REAL cpp   ! Cp
     25      REAL kappa ! kappa=R/Cp
     26      REAL cotot
     27      REAL unsim ! = 1./iim
     28      REAL g ! (m/s2) gravity
     29      REAL omeg ! (rad/s) rotation rate of the planet
    1830      REAL dissip_factz,dissip_deltaz,dissip_zref
    1931      INTEGER iflag_top_bound
    2032      REAL tau_top_bound
     33      REAL daylen ! length of solar day, in 'standard' day length
     34      REAL year_day ! Number of standard days in a year
     35      REAL molmass ! (g/mol) molar mass of the atmosphere
    2136
    2237

  • LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/friction.F

    r1403 r1437  
    66
    77      USE control_mod
    8  
     8#ifdef CPP_IOIPSL
     9      USE IOIPSL
     10#else
     11! if not using IOIPSL, we still need to use (a local version of) getin
     12      USE ioipsl_getincom
     13#endif
     14     
    915      IMPLICIT NONE
    1016
    11 c=======================================================================
    12 c
    13 c
    14 c   Objet:
    15 c   ------
    16 c
    17 c  ***********
    18 c    Friction
    19 c  ***********
    20 c
    21 c=======================================================================
     17!=======================================================================
     18!
     19!   Friction for the Newtonian case:
     20!   --------------------------------
     21!    2 possibilities (depending on flag 'friction_type'
     22!     friction_type=0 : A friction that is only applied to the lowermost
     23!                       atmospheric layer
     24!     friction_type=1 : Friction applied on all atmospheric layer (but
     25!       (default)       with stronger magnitude near the surface; see
     26!                       iniacademic.F)
     27!=======================================================================
    2228
    2329#include "dimensions.h"
    … …  
    2531#include "comgeom2.h"
    2632#include "comconst.h"
     33#include "iniprint.h"
     34#include "academic.h"
    2735
    28       REAL pdt
     36! arguments:
     37      REAL,INTENT(out) :: ucov( iip1,jjp1,llm )
     38      REAL,INTENT(out) :: vcov( iip1,jjm,llm )
     39      REAL,INTENT(in) :: pdt ! time step
     40
     41! local variables:
     42
    2943      REAL modv(iip1,jjp1),zco,zsi
    3044      REAL vpn,vps,upoln,upols,vpols,vpoln
    3145      REAL u2(iip1,jjp1),v2(iip1,jjm)
    32       REAL ucov( iip1,jjp1,llm ),vcov( iip1,jjm,llm )
    33       INTEGER  i,j
    34       REAL cfric
    35       parameter (cfric=1.e-5)
     46      INTEGER  i,j,l
     47      REAL,PARAMETER :: cfric=1.e-5
     48      LOGICAL,SAVE :: firstcall=.true.
     49      INTEGER,SAVE :: friction_type=1
     50      CHARACTER(len=20) :: modname="friction"
     51      CHARACTER(len=80) :: abort_message
     52     
     53      IF (firstcall) THEN
     54        ! set friction type
     55        call getin("friction_type",friction_type)
     56        if ((friction_type.lt.0).or.(friction_type.gt.1)) then
     57          abort_message="wrong friction type"
     58          write(lunout,*)'Friction: wrong friction type',friction_type
     59          call abort_gcm(modname,abort_message,42)
     60        endif
     61        firstcall=.false.
     62      ENDIF
    3663
    37 
     64      if (friction_type.eq.0) then
    3865c   calcul des composantes au carre du vent naturel
    3966      do j=1,jjp1
    … …  
    96123         vcov(iip1,j,1)=vcov(1,j,1)
    97124      enddo
     125      endif ! of if (friction_type.eq.0)
    98126
     127      if (friction_type.eq.1) then
     128        do l=1,llm
     129          ucov(:,:,l)=ucov(:,:,l)*(1.-pdt*kfrict(l))
     130          vcov(:,:,l)=vcov(:,:,l)*(1.-pdt*kfrict(l))
     131        enddo
     132      endif
     133     
    99134      RETURN
    100135      END

  • LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/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       

  • LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/leapfrog.F

    r1403 r1437  
    430430
    431431      IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
    432 c   Calcul academique de la physique = Rappel Newtonien + friction
    433 c   --------------------------------------------------------------
    434        teta(:,:)=teta(:,:)
    435      s  -iphysiq*dtvr*(teta(:,:)-tetarappel(:,:))/taurappel
    436        call friction(ucov,vcov,iphysiq*dtvr)
    437       ENDIF
     432!   Academic case : Simple friction and Newtonan relaxation
     433!   -------------------------------------------------------
     434        DO l=1,llm   
     435          DO ij=1,ip1jmp1
     436           teta(ij,l)=teta(ij,l)-dtvr*
     437     &      (teta(ij,l)-tetarappel(ij,l))*(knewt_g+knewt_t(l)*clat4(ij))
     438          ENDDO
     439        ENDDO ! of DO l=1,llm
     440          call friction(ucov,vcov,dtvr)
     441       
     442        ! Sponge layer (if any)
     443        IF (ok_strato) THEN
     444          dufi(:,:)=0.
     445          dvfi(:,:)=0.
     446          dtetafi(:,:)=0.
     447          dqfi(:,:,:)=0.
     448          dpfi(:)=0.
     449          CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
     450          CALL addfi( dtvr, leapf, forward   ,
     451     $                  ucov, vcov, teta , q   ,ps ,
     452     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
     453        ENDIF ! of IF (ok_strato)
     454      ENDIF ! of IF (iflag_phys.EQ.2)
    438455
    439456

  • LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/academic.h

    r774 r1437  
    11!
    2 ! $Header$
     2! $Id$
    33!
    4       real tetarappel(ip1jmp1,llm),taurappel
    5       common/academic/tetarappel,taurappel
     4      common/academic/tetarappel,knewt_t,kfrict,knewt_g,clat4
     5      real :: tetarappel(ip1jmp1,llm)
     6      real :: knewt_t(llm)
     7      real :: kfrict(llm)
     8      real :: knewt_g
     9      real :: clat4(ip1jmp1)

  • LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/comconst.h

    r1279 r1437  
    55! INCLUDE comconst.h
    66
    7       COMMON/comconst/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,           &
    8      & dtvr,daysec,                                                     &
     7      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
     8     &                 iflag_top_bound
     9      COMMON/comconstr/dtvr,daysec,                                     &
    910     & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg              &
    1011     &                   ,dissip_factz,dissip_deltaz,dissip_zref        &
    11      &                   ,iflag_top_bound,tau_top_bound
     12     &                   ,tau_top_bound,                                &
     13     & daylen,year_day,molmass
    1214
    1315
    1416      INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
    15       REAL dtvr,daysec
    16       REAL pi,dtphys,dtdiss,rad,r,cpp,kappa
    17       REAL cotot,unsim,g,omeg
     17      REAL dtvr ! dynamical time step (in s)
     18      REAL daysec !length (in s) of a standard day
     19      REAL pi    ! something like 3.14159....
     20      REAL dtphys ! (s) time step for the physics
     21      REAL dtdiss ! (s) time step for the dissipation
     22      REAL rad ! (m) radius of the planet
     23      REAL r ! Gas constant R=8.31 J.K-1.mol-1
     24      REAL cpp   ! Cp
     25      REAL kappa ! kappa=R/Cp
     26      REAL cotot
     27      REAL unsim ! = 1./iim
     28      REAL g ! (m/s2) gravity
     29      REAL omeg ! (rad/s) rotation rate of the planet
    1830      REAL dissip_factz,dissip_deltaz,dissip_zref
    1931      INTEGER iflag_top_bound
    2032      REAL tau_top_bound
     33      REAL daylen ! length of solar day, in 'standard' day length
     34      REAL year_day ! Number of standard days in a year
     35      REAL molmass ! (g/mol) molar mass of the atmosphere
    2136
    2237

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

    r1403 r1437  
    66      USE parallel
    77      USE control_mod
     8#ifdef CPP_IOIPSL
     9      USE IOIPSL
     10#else
     11! if not using IOIPSL, we still need to use (a local version of) getin
     12      USE ioipsl_getincom
     13#endif
    814      IMPLICIT NONE
    915
    10 c=======================================================================
    11 c
    12 c
    13 c   Objet:
    14 c   ------
    15 c
    16 c  ***********
    17 c    Friction
    18 c  ***********
    19 c
    20 c=======================================================================
     16!=======================================================================
     17!
     18!   Friction for the Newtonian case:
     19!   --------------------------------
     20!    2 possibilities (depending on flag 'friction_type'
     21!     friction_type=0 : A friction that is only applied to the lowermost
     22!                       atmospheric layer
     23!     friction_type=1 : Friction applied on all atmospheric layer (but
     24!       (default)       with stronger magnitude near the surface; see
     25!                       iniacademic.F)
     26!=======================================================================
    2127
    2228#include "dimensions.h"
    … …  
    2430#include "comgeom2.h"
    2531#include "comconst.h"
    26 
    27       REAL pdt
     32#include "iniprint.h"
     33#include "academic.h"
     34
     35! arguments:
     36      REAL,INTENT(out) :: ucov( iip1,jjp1,llm )
     37      REAL,INTENT(out) :: vcov( iip1,jjm,llm )
     38      REAL,INTENT(in) :: pdt ! time step
     39
     40! local variables:
    2841      REAL modv(iip1,jjp1),zco,zsi
    2942      REAL vpn,vps,upoln,upols,vpols,vpoln
    3043      REAL u2(iip1,jjp1),v2(iip1,jjm)
    31       REAL ucov( iip1,jjp1,llm ),vcov( iip1,jjm,llm )
    32       INTEGER  i,j
    33       REAL cfric
    34       parameter (cfric=1.e-5)
     44      INTEGER  i,j,l
     45      REAL,PARAMETER :: cfric=1.e-5
     46      LOGICAL,SAVE :: firstcall=.true.
     47      INTEGER,SAVE :: friction_type=1
     48      CHARACTER(len=20) :: modname="friction_p"
     49      CHARACTER(len=80) :: abort_message
     50!$OMP THREADPRIVATE(firstcall,friction_type)
    3551      integer :: jjb,jje
    3652
    37 
     53!$OMP SINGLE
     54      IF (firstcall) THEN
     55        ! set friction type
     56        call getin("friction_type",friction_type)
     57        if ((friction_type.lt.0).or.(friction_type.gt.1)) then
     58          abort_message="wrong friction type"
     59          write(lunout,*)'Friction: wrong friction type',friction_type
     60          call abort_gcm(modname,abort_message,42)
     61        endif
     62        firstcall=.false.
     63      ENDIF
     64!$OMP END SINGLE COPYPRIVATE(friction_type,firstcall)
     65
     66      if (friction_type.eq.0) then ! friction on first layer only
     67!$OMP SINGLE
    3868c   calcul des composantes au carre du vent naturel
    3969      jjb=jj_begin
    … …  
    138168         vcov(iip1,j,1)=vcov(1,j,1)
    139169      enddo
     170!$OMP END SINGLE
     171      endif ! of if (friction_type.eq.0)
     172
     173      if (friction_type.eq.1) then
     174       ! for ucov()
     175        jjb=jj_begin
     176        jje=jj_end
     177        if (pole_nord) jjb=jj_begin+1
     178        if (pole_sud) jje=jj_end-1
     179
     180!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     181        do l=1,llm
     182          ucov(1:iip1,jjb:jje,l)=ucov(1:iip1,jjb:jje,l)*
     183     &                                  (1.-pdt*kfrict(l))
     184        enddo
     185!$OMP END DO NOWAIT
     186       
     187       ! for vcoc()
     188        jjb=jj_begin
     189        jje=jj_end
     190        if (pole_sud) jje=jj_end-1
     191       
     192!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     193        do l=1,llm
     194          vcov(1:iip1,jjb:jje,l)=vcov(1:iip1,jjb:jje,l)*
     195     &                                  (1.-pdt*kfrict(l))
     196        enddo
     197!$OMP END DO
     198      endif ! of if (friction_type.eq.1)
    140199
    141200      RETURN

  • 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       

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

    r1403 r1437  
    968968
    969969      IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
    970 c   Calcul academique de la physique = Rappel Newtonien + fritcion
    971 c   --------------------------------------------------------------
    972 cym       teta(:,:)=teta(:,:)
    973 cym     s  -iphysiq*dtvr*(teta(:,:)-tetarappel(:,:))/taurappel
     970!   Academic case : Simple friction and Newtonan relaxation
     971!   -------------------------------------------------------
    974972       ijb=ij_begin
    975973       ije=ij_end
    976974!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    977975       do l=1,llm
    978        teta(ijb:ije,l)=teta(ijb:ije,l)
    979      &  -iphysiq*dtvr*(teta(ijb:ije,l)-tetarappel(ijb:ije,l))/taurappel
    980        enddo
     976        teta(ijb:ije,l)=teta(ijb:ije,l)-dtvr*
     977     &         (teta(ijb:ije,l)-tetarappel(ijb:ije,l))*
     978     &                  (knewt_g+knewt_t(l)*clat4(ijb:ije))
     979       enddo ! of do l=1,llm
    981980!$OMP END DO
    982981
    … …  
    987986       call WaitRequest(Request_Physic)     
    988987c$OMP BARRIER
    989 !$OMP MASTER
    990        call friction_p(ucov,vcov,iphysiq*dtvr)
    991 !$OMP END MASTER
     988       call friction_p(ucov,vcov,dtvr)
    992989!$OMP BARRIER
     990
     991        ! Sponge layer (if any)
     992        IF (ok_strato) THEN
     993          ! set dufi,dvfi,... to zero
     994          ijb=ij_begin
     995          ije=ij_end
     996!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     997          do l=1,llm
     998            dufi(ijb:ije,l)=0
     999            dtetafi(ijb:ije,l)=0
     1000            dqfi(ijb:ije,l,1:nqtot)=0
     1001          enddo
     1002!$OMP END DO
     1003!$OMP SINGLE
     1004          dpfi(ijb:ije)=0
     1005!$OMP END SINGLE
     1006          ijb=ij_begin
     1007          ije=ij_end
     1008          if (pole_sud) ije=ije-iip1
     1009!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     1010          do l=1,llm
     1011            dvfi(ijb:ije,l)=0
     1012          enddo
     1013!$OMP END DO
     1014
     1015          CALL top_bound_p(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
     1016          CALL addfi_p( dtvr, leapf, forward   ,
     1017     $                  ucov, vcov, teta , q   ,ps ,
     1018     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
     1019!$OMP BARRIER
     1020        ENDIF ! of IF (ok_strato)
    9931021      ENDIF ! of IF(iflag_phys.EQ.2)
    9941022
Note: See TracChangeset for help on using the changeset viewer.