Ignore:
Timestamp:
Nov 18, 2010, 1:01:24 PM (14 years ago)
Author:
Laurent Fairhead
Message:

Merge of LMDZ5V1.0-dev branch r1453 into LMDZ5 trunk r1434

Fusion entre la version r1453 de la branche de développement LMDZ5V1.0-dev
et le tronc LMDZ5 (r1434)

Location:
LMDZ5/trunk
Files:
15 edited
1 copied

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk

  • LMDZ5/trunk/libf/dyn3d/academic.h

    r524 r1454  
    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/trunk/libf/dyn3d/addfi.F

    r1146 r1454  
    11!
    2 ! $Header$
     2! $Id$
    33!
    44      SUBROUTINE addfi(pdt, leapf, forward,
    … …  
    77
    88      USE infotrac, ONLY : nqtot
     9      USE control_mod, ONLY : planet_type
    910      IMPLICIT NONE
    1011c
    … …  
    116117      ENDDO
    117118 
    118       DO iq = 1, 2
     119      if (planet_type=="earth") then
     120      ! earth case, special treatment for first 2 tracers (water)
     121       DO iq = 1, 2
    119122         DO k = 1,llm
    120123            DO j = 1,ip1jmp1
    … …  
    123126            ENDDO
    124127         ENDDO
    125       ENDDO
     128       ENDDO
    126129
    127       DO iq = 3, nqtot
     130       DO iq = 3, nqtot
    128131         DO k = 1,llm
    129132            DO j = 1,ip1jmp1
    … …  
    132135            ENDDO
    133136         ENDDO
    134       ENDDO
     137       ENDDO
     138      else
     139      ! general case, treat all tracers equally)
     140       DO iq = 1, nqtot
     141         DO k = 1,llm
     142            DO j = 1,ip1jmp1
     143               pq(j,k,iq)= pq(j,k,iq) + pdqfi(j,k,iq) * pdt
     144               pq(j,k,iq)= AMAX1( pq(j,k,iq), qtestt )
     145            ENDDO
     146         ENDDO
     147       ENDDO
     148      endif ! of if (planet_type=="earth")
    135149
    136150

  • LMDZ5/trunk/libf/dyn3d/advtrac.F

    r1403 r1454  
    236236            call vlsplt(q(1,1,iq),2.,massem,wg,pbarug,pbarvg,dtvr)
    237237
    238 
    239238c   ----------------------------------------------------------------
    240239c   Schema "pseudo amont" + test sur humidite specifique

  • LMDZ5/trunk/libf/dyn3d/caladvtrac.F

    r1403 r1454  
    88     *                   flxw, pk)
    99c
    10       USE infotrac
    11       USE control_mod
     10      USE infotrac, ONLY : nqtot
     11      USE control_mod, ONLY : iapp_tracvl,planet_type
    1212 
    1313      IMPLICIT NONE
    … …  
    3030c   ----------
    3131      REAL pbaru( ip1jmp1,llm ),pbarv( ip1jm,llm),masse(ip1jmp1,llm)
    32       REAL p( ip1jmp1,llmp1),q( ip1jmp1,llm,nqtot),dq( ip1jmp1,llm,2 )
     32      REAL p( ip1jmp1,llmp1),q( ip1jmp1,llm,nqtot)
     33      real :: dq(ip1jmp1,llm,nqtot)
    3334      REAL teta( ip1jmp1,llm),pk( ip1jmp1,llm)
    3435      REAL               :: flxw(ip1jmp1,llm)
    … …  
    4950cc
    5051c
     52! Earth-specific stuff for the first 2 tracers (water)
     53      if (planet_type.eq."earth") then
    5154C initialisation
    52         dq = 0.
    53 
    54         CALL SCOPY( 2 * ijp1llm, q, 1, dq, 1 )
    55 
     55        dq(:,:,1:2)=q(:,:,1:2)
     56       
    5657c  test des valeurs minmax
    5758cc        CALL minmaxq(q(1,1,1),1.e33,-1.e33,'Eau vapeur (a) ')
    5859cc        CALL minmaxq(q(1,1,2),1.e33,-1.e33,'Eau liquide(a) ')
    59 
     60      endif ! of if (planet_type.eq."earth")
    6061c   advection
    6162
    … …  
    6364     *       p,  masse,q,iapptrac, teta,
    6465     .       flxw, pk)
     66
    6567c
    6668
    67          IF( iapptrac.EQ.iapp_tracvl ) THEN
     69      IF( iapptrac.EQ.iapp_tracvl ) THEN
     70        if (planet_type.eq."earth") then
     71! Earth-specific treatment for the first 2 tracers (water)
    6872c
    6973cc          CALL minmaxq(q(1,1,1),1.e33,-1.e33,'Eau vapeur     ')
    … …  
    7882          ENDDO
    7983         
    80           if (planet_type.eq."earth") then
    81 ! Earth-specific treatment of first 2 tracers (water)
    82             CALL qminimum( q, 2, finmasse )
    83           endif
     84          CALL qminimum( q, 2, finmasse )
    8485
    8586          CALL SCOPY   ( ip1jmp1*llm, masse, 1, finmasse,       1 )
    … …  
    100101           ENDDO
    101102c
    102          ELSE
    103            DO iq = 1 , 2
    104            DO l  = 1, llm
    105              DO ij = 1,ip1jmp1
    106               dq(ij,l,iq)  = 0.
    107              ENDDO
    108            ENDDO
    109            ENDDO
     103        endif ! of if (planet_type.eq."earth")
     104      ELSE
     105        if (planet_type.eq."earth") then
     106! Earth-specific treatment for the first 2 tracers (water)
     107          dq(:,:,1:2)=0.
     108        endif ! of if (planet_type.eq."earth")
     109      ENDIF ! of IF( iapptrac.EQ.iapp_tracvl )
    110110
    111 
    112          ENDIF
    113 
    114 c
    115 
    116 c  ... On appelle  qminimum uniquement  pour l'eau vapeur et liquide  ..
    117 
    118  
    119       RETURN
    120111      END
    121112

  • LMDZ5/trunk/libf/dyn3d/comconst.h

    r1279 r1454  
    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/trunk/libf/dyn3d/ener.h

    r524 r1454  
    11!
    2 ! $Header$
     2! $Id$
    33!
    4 c-----------------------------------------------------------------------
    5 c INCLUDE 'ener.h'
     4!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
     5!                 veillez à n'utiliser que des ! pour les commentaires
     6!                 et à bien positionner les & des lignes de continuation
     7!                 (les placer en colonne 6 et en colonne 73)
     8!
     9! INCLUDE 'ener.h'
    610
    7       COMMON/ener/ang0,etot0,ptot0,ztot0,stot0,
    8      *            ang,etot,ptot,ztot,stot,rmsdpdt ,
    9      *            rmsv,gtot(llmm1)
     11      COMMON/ener/ang0,etot0,ptot0,ztot0,stot0,                         &
     12     &            ang,etot,ptot,ztot,stot,rmsdpdt ,                     &
     13     &            rmsv,gtot(llmm1)
    1014
    11       REAL ang0,etot0,ptot0,ztot0,stot0,
    12      s     ang,etot,ptot,ztot,stot,rmsdpdt,rmsv,gtot
     15      REAL ang0,etot0,ptot0,ztot0,stot0,                                &
     16     &     ang,etot,ptot,ztot,stot,rmsdpdt,rmsv,gtot
    1317
    14 c-----------------------------------------------------------------------
     18!-----------------------------------------------------------------------

  • LMDZ5/trunk/libf/dyn3d/friction.F

    r1403 r1454  
    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/trunk/libf/dyn3d/gcm.F

    r1403 r1454  
    249249        endif
    250250
    251         if (planet_type.eq."earth") then
    252 #ifdef CPP_EARTH
     251!        if (planet_type.eq."earth") then
    253252! Load an Earth-format start file
    254253         CALL dynetat0("start.nc",vcov,ucov,
    255254     &              teta,q,masse,ps,phis, time_0)
    256 #else
    257         ! SW model also has Earth-format start files
    258         ! (but can be used without the CPP_EARTH directive)
    259           if (iflag_phys.eq.0) then
    260             CALL dynetat0("start.nc",vcov,ucov,
    261      &              teta,q,masse,ps,phis, time_0)
    262           endif
    263 #endif
    264         endif ! of if (planet_type.eq."earth")
     255!        endif ! of if (planet_type.eq."earth")
    265256       
    266257c       write(73,*) 'ucov',ucov
    … …  
    468459#endif
    469460
    470       if (planet_type.eq."earth") then
     461!      if (planet_type.eq."earth") then
     462! Write an Earth-format restart file
    471463        CALL dynredem0("restart.nc", day_end, phis)
    472       endif
     464!      endif
    473465
    474466      ecripar = .TRUE.

  • LMDZ5/trunk/libf/dyn3d/grid_noro.F

    r1403 r1454  
    458458C MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
    459459
    460       PARAMETER (ISMo=400,JSMo=200)
    461       REAL X(IMAR,JMAR),XF(ISMo,JSMo)
     460      REAL X(IMAR,JMAR),XF(IMAR,JMAR)
    462461      real WEIGHTpb(-1:1,-1:1)
    463462
    464       if(imar.gt.ismo) stop'surdimensionner ismo dans mva9 (grid_noro)'
    465       if(jmar.gt.jsmo) stop'surdimensionner jsmo dans mva9 (grid_noro)'
    466      
     463
    467464      SUM=0.
    468465      DO IS=-1,1

  • LMDZ5/trunk/libf/dyn3d/infotrac.F90

    r1403 r1454  
    3131
    3232  SUBROUTINE infotrac_init
    33 
    3433    USE control_mod
    35  
    3634    IMPLICIT NONE
    3735!=======================================================================
    … …  
    6361    CHARACTER(len=1), DIMENSION(3)  :: txts
    6462    CHARACTER(len=2), DIMENSION(9)  :: txtp
    65     CHARACTER(len=13)               :: str1,str2
     63    CHARACTER(len=23)               :: str1,str2
    6664 
    6765    INTEGER :: nqtrue  ! number of tracers read from tracer.def, without higer order of moment
    6866    INTEGER :: iq, new_iq, iiq, jq, ierr
    69     INTEGER, EXTERNAL :: lnblnk
    70  
     67
     68    character(len=*),parameter :: modname="infotrac_init"
    7169!-----------------------------------------------------------------------
    7270! Initialization :
    … …  
    102100       OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
    103101       IF(ierr.EQ.0) THEN
    104           WRITE(lunout,*) 'Open traceur.def : ok'
     102          WRITE(lunout,*) trim(modname),': Open traceur.def : ok'
    105103          READ(90,*) nqtrue
    106104       ELSE
    107           WRITE(lunout,*) 'Problem in opening traceur.def'
    108           WRITE(lunout,*) 'ATTENTION using defaut values'
    109           nqtrue=4 ! Defaut value
    110        END IF
    111        ! Attention! Only for planet_type=='earth'
    112        nbtr=nqtrue-2
     105          WRITE(lunout,*) trim(modname),': Problem in opening traceur.def'
     106          WRITE(lunout,*) trim(modname),': WARNING using defaut values'
     107          if (planet_type=='earth') then
     108            nqtrue=4 ! Default value for Earth
     109          else
     110            nqtrue=1 ! Default value for other planets
     111          endif
     112       END IF
     113       if ( planet_type=='earth') then
     114         ! For Earth, water vapour & liquid tracers are not in the physics
     115         nbtr=nqtrue-2
     116       else
     117         ! Other planets (for now); we have the same number of tracers
     118         ! in the dynamics than in the physics
     119         nbtr=nqtrue
     120       endif
    113121    ELSE
    114122       ! nbtr has been read from INCA by init_cont_lmdz() in gcm.F
    … …  
    116124    END IF
    117125
    118     IF (nqtrue < 2) THEN
    119        WRITE(lunout,*) 'nqtrue=',nqtrue, ' is not allowded. 2 tracers is the minimum'
     126    IF ((planet_type=="earth").and.(nqtrue < 2)) THEN
     127       WRITE(lunout,*) trim(modname),': nqtrue=',nqtrue, ' is not allowded. 2 tracers is the minimum'
    120128       CALL abort_gcm('infotrac_init','Not enough tracers',1)
    121129    END IF
    … …  
    158166          ! Continue to read tracer.def
    159167          DO iq=1,nqtrue
    160              READ(90,999) hadv(iq),vadv(iq),tnom_0(iq)
     168             READ(90,*) hadv(iq),vadv(iq),tnom_0(iq)
    161169          END DO
    162170          CLOSE(90) 
    163        ELSE ! Without tracer.def
     171       ELSE ! Without tracer.def, set default values
     172         if (planet_type=="earth") then
     173          ! for Earth, default is to have 4 tracers
    164174          hadv(1) = 14
    165175          vadv(1) = 14
    … …  
    174184          vadv(4) = 10
    175185          tnom_0(4) = 'PB'
     186         else ! default for other planets
     187          hadv(1) = 10
     188          vadv(1) = 10
     189          tnom_0(1) = 'dummy'
     190         endif ! of if (planet_type=="earth")
    176191       END IF
    177192       
    178        WRITE(lunout,*) 'Valeur de traceur.def :'
    179        WRITE(lunout,*) 'nombre de traceurs ',nqtrue
     193       WRITE(lunout,*) trim(modname),': Valeur de traceur.def :'
     194       WRITE(lunout,*) trim(modname),': nombre de traceurs ',nqtrue
    180195       DO iq=1,nqtrue
    181196          WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq)
    … …  
    219234          new_iq=new_iq+10 ! 9 tracers added
    220235       ELSE
    221           WRITE(lunout,*) 'This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
     236          WRITE(lunout,*) trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
    222237          CALL abort_gcm('infotrac_init','Bad choice of advection schema - 1',1)
    223238       END IF
    … …  
    229244       nqtot = new_iq
    230245
    231        WRITE(lunout,*) 'The choice of advection schema for one or more tracers'
     246       WRITE(lunout,*) trim(modname),': The choice of advection schema for one or more tracers'
    232247       WRITE(lunout,*) 'makes it necessary to add tracers'
    233        WRITE(lunout,*) nqtrue,' is the number of true tracers'
    234        WRITE(lunout,*) nqtot, ' is the total number of tracers needed'
     248       WRITE(lunout,*) trim(modname)//': ',nqtrue,' is the number of true tracers'
     249       WRITE(lunout,*) trim(modname)//': ',nqtot, ' is the total number of tracers needed'
    235250
    236251    ELSE
    … …  
    260275          iadv(new_iq)=11
    261276       ELSE
    262           WRITE(lunout,*)'This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
     277          WRITE(lunout,*)trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
     278
    263279          CALL abort_gcm('infotrac_init','Bad choice of advection schema - 2',1)
    264280       END IF
    … …  
    267283       tname(new_iq)= tnom_0(iq)
    268284       IF (iadv(new_iq)==0) THEN
    269           ttext(new_iq)=str1(1:lnblnk(str1))
     285          ttext(new_iq)=trim(str1)
    270286       ELSE
    271           ttext(new_iq)=str1(1:lnblnk(str1))//descrq(iadv(new_iq))
     287          ttext(new_iq)=trim(tnom_0(iq))//descrq(iadv(new_iq))
    272288       END IF
    273289
    … …  
    278294             new_iq=new_iq+1
    279295             iadv(new_iq)=-20
    280              ttext(new_iq)=str2(1:lnblnk(str2))//txts(jq)
    281              tname(new_iq)=str1(1:lnblnk(str1))//txts(jq)
     296             ttext(new_iq)=trim(str2)//txts(jq)
     297             tname(new_iq)=trim(str1)//txts(jq)
    282298          END DO
    283299       ELSE IF (iadv(new_iq)==30) THEN
    … …  
    285301             new_iq=new_iq+1
    286302             iadv(new_iq)=-30
    287              ttext(new_iq)=str2(1:lnblnk(str2))//txtp(jq)
    288              tname(new_iq)=str1(1:lnblnk(str1))//txtp(jq)
     303             ttext(new_iq)=trim(str2)//txtp(jq)
     304             tname(new_iq)=trim(str1)//txtp(jq)
    289305          END DO
    290306       END IF
    … …  
    305321
    306322
    307     WRITE(lunout,*) 'Information stored in infotrac :'
    308     WRITE(lunout,*) 'iadv  niadv tname  ttext :'
     323    WRITE(lunout,*) trim(modname),': Information stored in infotrac :'
     324    WRITE(lunout,*) trim(modname),': iadv  niadv tname  ttext :'
    309325    DO iq=1,nqtot
    310        WRITE(lunout,*) iadv(iq),niadv(iq), tname(iq), ttext(iq)
     326       WRITE(lunout,*) iadv(iq),niadv(iq),&
     327       ' ',trim(tname(iq)),' ',trim(ttext(iq))
    311328    END DO
    312329
    … …  
    317334    DO iq=1,nqtot
    318335       IF (iadv(iq)/=10 .AND. iadv(iq)/=14 .AND. iadv(iq)/=0) THEN
    319           WRITE(lunout,*)'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
     336          WRITE(lunout,*)trim(modname),' STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
    320337          CALL abort_gcm('infotrac_init','In this version only iadv=10 and iadv=14 is tested!',1)
    321338       ELSE IF (iadv(iq)==14 .AND. iq/=1) THEN
    322           WRITE(lunout,*)'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
     339          WRITE(lunout,*)trim(modname),'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ'
    323340          CALL abort_gcm('infotrac_init','In this version iadv=14 is only permitted for water vapour!',1)
    324341       END IF
    … …  
    331348    DEALLOCATE(tracnam)
    332349
    333 999 FORMAT (i2,1x,i2,1x,a15)
    334 
    335350  END SUBROUTINE infotrac_init
    336351

  • LMDZ5/trunk/libf/dyn3d/iniacademic.F

    r1403 r1454  
    88      USE filtreg_mod
    99      USE infotrac, ONLY : nqtot
    10       USE control_mod
     10      USE control_mod, ONLY: day_step,planet_type
     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
    … …  
    8295! 1. Initializations for Earth-like case
    8396! --------------------------------------
    84       if (planet_type=="earth") then
    85 c
     97c
     98        ! initialize planet radius, rotation rate,...
     99        call conf_planete
     100
    86101        time_0=0.
    87         day_ref=0
     102        day_ref=1
    88103        annee_ref=0
    89104
    90105        im         = iim
    91106        jm         = jjm
    92         day_ini    = 0
    93         omeg       = 4.*asin(1.)/86400.
    94         rad    = 6371229.
    95         g      = 9.8
    96         daysec = 86400.
     107        day_ini    = 1
    97108        dtvr    = daysec/REAL(day_step)
    98109        zdtvr=dtvr
    99         kappa  = 0.2857143
    100         cpp    = 1004.70885
    101         preff     = 101325.
    102         pa        =  50000.
    103110        etot0      = 0.
    104111        ptot0      = 0.
    … …  
    120127          if (.not.read_start) then
    121128            phis(:)=0.
    122             q(:,:,1)=1.e-10
    123             q(:,:,2)=1.e-15
    124             q(:,:,3:nqtot)=0.
     129            q(:,:,:)=0
    125130            CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
    126131          endif
    … …  
    129134        if (iflag_phys.eq.2) then
    130135          ! 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 
     136         
     137!         if (planet_type=="earth") then
     138
     139          ! 1. local parameters
     140          ! by convention, winter is in the southern hemisphere
     141          ! Geostrophic wind or no wind?
     142          ok_geost=.TRUE.
     143          CALL getin('ok_geost',ok_geost)
     144          ! Constants for Newtonian relaxation and friction
     145          k_f=1.                !friction
     146          CALL getin('k_j',k_f)
     147          k_f=1./(daysec*k_f)
     148          k_c_s=4.  !cooling surface
     149          CALL getin('k_c_s',k_c_s)
     150          k_c_s=1./(daysec*k_c_s)
     151          k_c_a=40. !cooling free atm
     152          CALL getin('k_c_a',k_c_a)
     153          k_c_a=1./(daysec*k_c_a)
     154          ! Constants for Teta equilibrium profile
     155          teta0=315.     ! mean Teta (S.H. 315K)
     156          CALL getin('teta0',teta0)
     157          ttp=200.       ! Tropopause temperature (S.H. 200K)
     158          CALL getin('ttp',ttp)
     159          eps=0.         ! Deviation to N-S symmetry(~0-20K)
     160          CALL getin('eps',eps)
     161          delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
     162          CALL getin('delt_y',delt_y)
     163          delt_z=10.     ! Vertical Gradient (S.H. 10K)
     164          CALL getin('delt_z',delt_z)
     165          ! Polar vortex
     166          ok_pv=.false.
     167          CALL getin('ok_pv',ok_pv)
     168          phi_pv=-50.            ! Latitude of edge of vortex
     169          CALL getin('phi_pv',phi_pv)
     170          phi_pv=phi_pv*pi/180.
     171          dphi_pv=5.             ! Width of the edge
     172          CALL getin('dphi_pv',dphi_pv)
     173          dphi_pv=dphi_pv*pi/180.
     174          gam_pv=4.              ! -dT/dz vortex (in K/km)
     175          CALL getin('gam_pv',gam_pv)
     176         
     177          ! 2. Initialize fields towards which to relax
     178          ! Friction
     179          knewt_g=k_c_a
    141180          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)
     181           zsig=presnivs(l)/preff
     182           knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
     183           kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
     184          ENDDO
     185          DO j=1,jjp1
     186            clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
     187          ENDDO
     188         
     189          ! Potential temperature
     190          DO l=1,llm
     191           zsig=presnivs(l)/preff
     192           tetastrat=ttp*zsig**(-kappa)
     193           tetapv=tetastrat
     194           IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
     195               tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
     196           ENDIF
     197           DO j=1,jjp1
     198             ! Troposphere
     199             ddsin=sin(rlatu(j))
     200             tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin
     201     &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
     202             ! Profil stratospherique isotherme (+vortex)
     203             w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
     204             tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
     205             tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
     206           ENDDO
    156207          ENDDO ! of DO l=1,llm
     208 
     209!          CALL writefield('theta_eq',tetajl)
    157210
    158211          do l=1,llm
    … …  
    165218          enddo
    166219
    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))
     220
     221!         else
     222!          write(lunout,*)"iniacademic: planet types other than earth",
     223!     &                   " not implemented (yet)."
     224!          stop
     225!         endif ! of if (planet_type=="earth")
     226
     227          ! 3. Initialize fields (if necessary)
     228          IF (.NOT. read_start) THEN
     229            ! surface pressure
     230            ps(:)=preff
     231            ! ground geopotential
     232            phis(:)=0.
     233           
     234            CALL pression ( ip1jmp1, ap, bp, ps, p       )
     235            CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     236            CALL massdair(p,masse)
     237
     238            ! bulk initialization of temperature
     239            teta(:,:)=tetarappel(:,:)
     240           
     241            ! geopotential
     242            CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
     243           
     244            ! winds
     245            if (ok_geost) then
     246              call ugeostr(phi,ucov)
     247            else
     248              ucov(:,:)=0.
     249            endif
     250            vcov(:,:)=0.
     251           
     252            ! bulk initialization of tracers
     253            if (planet_type=="earth") then
     254              ! Earth: first two tracers will be water
     255              do i=1,nqtot
     256                if (i.eq.1) q(:,:,i)=1.e-10
     257                if (i.eq.2) q(:,:,i)=1.e-15
     258                if (i.gt.2) q(:,:,i)=0.
     259              enddo
     260            else
     261              q(:,:,:)=0
     262            endif ! of if (planet_type=="earth")
     263
     264            ! add random perturbation to temperature
     265            idum  = -1
     266            zz = ran1(idum)
     267            idum  = 0
     268            do l=1,llm
     269              do ij=iip2,ip1jm
     270                teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
     271              enddo
    190272            enddo
    191           enddo
    192 
    193           do l=1,llm
    194             do ij=1,ip1jmp1,iip1
    195               teta(ij+iim,l)=teta(ij,l)
     273
     274            ! maintain periodicity in longitude
     275            do l=1,llm
     276              do ij=1,ip1jmp1,iip1
     277                teta(ij+iim,l)=teta(ij,l)
     278              enddo
    196279            enddo
    197           enddo
    198 
    199 
    200 
    201 c     PRINT *,' Appel test_period avec tetarappel '
    202 c     CALL  test_period ( ucov,vcov,tetarappel,q,p,phis )
    203 c     PRINT *,' Appel test_period avec teta '
    204 c     CALL  test_period ( ucov,vcov,teta,q,p,phis )
    205 
    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.
     280
     281          ENDIF ! of IF (.NOT. read_start)
    211282        endif ! of if (iflag_phys.eq.2)
    212283       
    213       else
    214         write(lunout,*)"iniacademic: planet types other than earth",
    215      &                 " not implemented (yet)."
    216         stop
    217       endif ! of if (planet_type=="earth")
    218       return
    219284      END
    220285c-----------------------------------------------------------------------

  • LMDZ5/trunk/libf/dyn3d/integrd.F

    r1403 r1454  
    66     $     dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis,finvmaold )
    77
    8       USE control_mod
     8      use control_mod, only : planet_type
    99
    1010      IMPLICIT NONE
    … …  
    8181      CALL SCOPY(ip1jmp1*llm, masse, 1, massescr, 1)
    8282
    83       DO 2 ij = 1,ip1jmp1
     83      DO ij = 1,ip1jmp1
    8484       pscr (ij)    = ps(ij)
    8585       ps (ij)      = psm1(ij) + dt * dp(ij)
    86    2  CONTINUE
     86      ENDDO
    8787c
    8888      DO ij = 1,ip1jmp1
    … …  
    115115c    ............   integration  de  ucov, vcov,  h     ..............
    116116
    117       DO 10 l = 1,llm
    118 
    119       DO 4 ij = iip2,ip1jm
    120       uscr( ij )   =  ucov( ij,l )
    121       ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
    122    4  CONTINUE
    123 
    124       DO 5 ij = 1,ip1jm
    125       vscr( ij )   =  vcov( ij,l )
    126       vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
    127    5  CONTINUE
    128 
    129       DO 6 ij = 1,ip1jmp1
    130       hscr( ij )    =  teta(ij,l)
    131       teta ( ij,l ) = tetam1(ij,l) *  massem1(ij,l) / masse(ij,l)
    132      $                + dt * dteta(ij,l) / masse(ij,l)
    133    6  CONTINUE
     117      DO l = 1,llm
     118
     119       DO ij = iip2,ip1jm
     120        uscr( ij )   =  ucov( ij,l )
     121        ucov( ij,l ) = ucovm1( ij,l ) + dt * du( ij,l )
     122       ENDDO
     123
     124       DO ij = 1,ip1jm
     125        vscr( ij )   =  vcov( ij,l )
     126        vcov( ij,l ) = vcovm1( ij,l ) + dt * dv( ij,l )
     127       ENDDO
     128
     129       DO ij = 1,ip1jmp1
     130        hscr( ij )    =  teta(ij,l)
     131        teta ( ij,l ) = tetam1(ij,l) *  massem1(ij,l) / masse(ij,l)
     132     &                + dt * dteta(ij,l) / masse(ij,l)
     133       ENDDO
    134134
    135135c   ....  Calcul de la valeur moyenne, unique  aux poles pour  teta    ......
    136136c
    137137c
    138       DO  ij   = 1, iim
     138       DO  ij   = 1, iim
    139139        tppn(ij) = aire(   ij   ) * teta(  ij    ,l)
    140140        tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
    141       ENDDO
     141       ENDDO
    142142        tpn      = SSUM(iim,tppn,1)/apoln
    143143        tps      = SSUM(iim,tpps,1)/apols
    144144
    145       DO ij   = 1, iip1
     145       DO ij   = 1, iip1
    146146        teta(   ij   ,l)  = tpn
    147147        teta(ij+ip1jm,l)  = tps
    148       ENDDO
    149 c
    150 
    151       IF(leapf)  THEN
     148       ENDDO
     149c
     150
     151       IF(leapf)  THEN
    152152         CALL SCOPY ( ip1jmp1, uscr(1), 1, ucovm1(1, l), 1 )
    153153         CALL SCOPY (   ip1jm, vscr(1), 1, vcovm1(1, l), 1 )
    154154         CALL SCOPY ( ip1jmp1, hscr(1), 1, tetam1(1, l), 1 )
    155       END IF
    156 
    157   10  CONTINUE
     155       END IF
     156
     157      ENDDO ! of DO l = 1,llm
    158158
    159159
    … …  
    185185c$$$      ENDIF
    186186
    187         if (planet_type.eq."earth") then
     187      if (planet_type.eq."earth") then
    188188! Earth-specific treatment of first 2 tracers (water)
    189           DO l = 1, llm
    190            DO ij = 1, ip1jmp1
     189        DO l = 1, llm
     190          DO ij = 1, ip1jmp1
    191191            deltap(ij,l) =  p(ij,l) - p(ij,l+1)
    192            ENDDO
    193192          ENDDO
    194 
    195           CALL qminimum( q, nq, deltap )
    196          endif ! of if (planet_type.eq."earth")
     193        ENDDO
     194
     195        CALL qminimum( q, nq, deltap )
    197196
    198197c
    … …  
    200199c
    201200
    202       DO iq = 1, nq
     201       DO iq = 1, nq
    203202        DO l = 1, llm
    204203
    … …  
    216215
    217216        ENDDO
    218       ENDDO
    219 
    220 
    221          CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
     217       ENDDO
     218
     219
     220      CALL  SCOPY( ijp1llm , finvmasse, 1, finvmaold, 1 )
     221
     222      endif ! of if (planet_type.eq."earth")
    222223c
    223224c
    224225c     .....   FIN  de l'integration  de   q    .......
    225 
    226 15    continue
    227226
    228227c    .................................................................

  • LMDZ5/trunk/libf/dyn3d/leapfrog.F

    r1403 r1454  
    209209c   --------------------------------------------------
    210210
    211       dq=0.
     211      dq(:,:,:)=0.
    212212      CALL pression ( ip1jmp1, ap, bp, ps, p       )
    213213      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    … …  
    253253      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
    254254      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
    255 
    256 ! Ehouarn: what is this for? zqmin & zqmax are not used anyway ...
    257 !      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
    258255
    259256   2  CONTINUE
    … …  
    430427
    431428      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
     429!   Academic case : Simple friction and Newtonan relaxation
     430!   -------------------------------------------------------
     431        DO l=1,llm   
     432          DO ij=1,ip1jmp1
     433           teta(ij,l)=teta(ij,l)-dtvr*
     434     &      (teta(ij,l)-tetarappel(ij,l))*(knewt_g+knewt_t(l)*clat4(ij))
     435          ENDDO
     436        ENDDO ! of DO l=1,llm
     437          call friction(ucov,vcov,dtvr)
     438       
     439        ! Sponge layer (if any)
     440        IF (ok_strato) THEN
     441          dufi(:,:)=0.
     442          dvfi(:,:)=0.
     443          dtetafi(:,:)=0.
     444          dqfi(:,:,:)=0.
     445          dpfi(:)=0.
     446          CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
     447          CALL addfi( dtvr, leapf, forward   ,
     448     $                  ucov, vcov, teta , q   ,ps ,
     449     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
     450        ENDIF ! of IF (ok_strato)
     451      ENDIF ! of IF (iflag_phys.EQ.2)
    438452
    439453
    … …  
    615629
    616630
    617               if (planet_type.eq."earth") then
     631!              if (planet_type.eq."earth") then
    618632! Write an Earth-format restart file
    619633                CALL dynredem1("restart.nc",0.0,
    620634     &                         vcov,ucov,teta,q,masse,ps)
    621               endif ! of if (planet_type.eq."earth")
     635!              endif ! of if (planet_type.eq."earth")
    622636
    623637              CLOSE(99)
    … …  
    727741
    728742              IF(itau.EQ.itaufin) THEN
    729                 if (planet_type.eq."earth") then
     743!                if (planet_type.eq."earth") then
    730744                  CALL dynredem1("restart.nc",0.0,
    731745     &                           vcov,ucov,teta,q,masse,ps)
    732                 endif ! of if (planet_type.eq."earth")
     746!                endif ! of if (planet_type.eq."earth")
    733747              ENDIF ! of IF(itau.EQ.itaufin)
    734748

  • LMDZ5/trunk/libf/dyn3d/limit_netcdf.F90

    r1425 r1454  
    9797  kappa = 0.2857143
    9898  cpp   = 1004.70885
    99   dtvr  = daysec/FLOAT(day_step)
     99  dtvr  = daysec/REAL(day_step)
    100100  CALL inigeom
    101101
    … …  
    265265
    266266  DEALLOCATE(pctsrf_t,phy_sst,phy_bil,phy_alb,phy_rug)
    267 #endif
    268 ! of #ifdef CPP_EARTH
    269267
    270268
    … …  
    592590
    593591!--- Mid-months times
    594   mid_months(1)=0.5*FLOAT(mnth(1))
     592  mid_months(1)=0.5*REAL(mnth(1))
    595593  DO k=2,nm
    596     mid_months(k)=mid_months(k-1)+0.5*FLOAT(mnth(k-1)+mnth(k))
     594    mid_months(k)=mid_months(k-1)+0.5*REAL(mnth(k-1)+mnth(k))
    597595  END DO
    598596
    … …  
    626624!-------------------------------------------------------------------------------
    627625
     626#endif
     627! of #ifdef CPP_EARTH
    628628
    629629END SUBROUTINE limit_netcdf
Note: See TracChangeset for help on using the changeset viewer.