Changeset 1454 for LMDZ5/trunk/libf


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:
2 deleted
45 edited
4 copied

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk

  • LMDZ5/trunk/libf/cosp/icarus.F

    r1414 r1454  
    450450              write (6,*) 'rangevec:'
    451451              write (6,*) rangevec
    452               call flush(6)
    453452              STOP
    454453        endif

  • 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

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

    r774 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/dyn3dpar/addfi_p.F

    r1146 r1454  
    11!
    2 ! $Header$
     2! $Id$
    33!
    44      SUBROUTINE addfi_p(pdt, leapf, forward,
    … …  
    77      USE parallel
    88      USE infotrac, ONLY : nqtot
     9      USE control_mod, ONLY : planet_type
    910      IMPLICIT NONE
    1011c
    … …  
    154155c$OMP END MASTER
    155156 
    156       DO iq = 1, 2
     157      if (planet_type=="earth") then
     158      ! earth case, special treatment for first 2 tracers (water)
     159       DO iq = 1, 2
    157160c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    158161         DO k = 1,llm
    … …  
    163166         ENDDO
    164167c$OMP END DO NOWAIT
    165       ENDDO
    166 
    167       DO iq = 3, nqtot
     168       ENDDO
     169
     170       DO iq = 3, nqtot
    168171c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    169172         DO k = 1,llm
    … …  
    174177         ENDDO
    175178c$OMP END DO NOWAIT
    176       ENDDO
     179       ENDDO
     180      else
     181      ! general case, treat all tracers equally)
     182       DO iq = 1, nqtot
     183c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     184         DO k = 1,llm
     185            DO j = ijb,ije
     186               pq(j,k,iq)= pq(j,k,iq) + pdqfi(j,k,iq) * pdt
     187               pq(j,k,iq)= AMAX1( pq(j,k,iq), qtestt )
     188            ENDDO
     189         ENDDO
     190c$OMP END DO NOWAIT
     191       ENDDO
     192      endif ! of if (planet_type=="earth")
    177193
    178194c$OMP MASTER

  • LMDZ5/trunk/libf/dyn3dpar/advtrac_p.F

    r1403 r1454  
    132132ccc         CALL filtreg ( massem ,jjp1, llm,-2, 2, .TRUE., 1 )
    133133c
    134       ENDIF
     134      ENDIF ! of IF(iadvtr.EQ.0)
    135135
    136136      iadvtr   = iadvtr+1
    … …  
    266266cym      ----> Revérifier lors de la parallélisation des autres schemas
    267267   
    268 cym          call massbar_p(massem,massebx,masseby)         
     268cym          call massbar_p(massem,massebx,masseby) 
    269269
    270270          call vlspltgen_p( q,iadv, 2., massem, wg ,
    … …  
    452452c$OMP BARRIER
    453453
    454       ijb=ij_begin
    455       ije=ij_end
     454      if (planet_type=="earth") then
     455
     456        ijb=ij_begin
     457        ije=ij_end
    456458
    457459c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    458        DO l = 1, llm
     460        DO l = 1, llm
    459461         DO ij = ijb, ije
    460462           finmasse(ij,l) =  p(ij,l) - p(ij,l+1)
    461463         ENDDO
    462        ENDDO
     464        ENDDO
    463465c$OMP END DO
    464466
    465        CALL qminimum_p( q, 2, finmasse )
     467        CALL qminimum_p( q, 2, finmasse )
    466468
    467469c------------------------------------------------------------------
    … …  
    496498c$OMP BARRIER   
    497499          iadvtr=0
     500      endif ! of if (planet_type=="earth")
    498501       ENDIF ! if iadvtr.EQ.iapp_tracvl
    499502

  • LMDZ5/trunk/libf/dyn3dpar/caladvtrac_p.F

    r1403 r1454  
    88     *                   flxw, pk, iapptrac)
    99      USE parallel
    10       USE infotrac
    11       USE control_mod
     10      USE infotrac, ONLY : nqtot
     11      USE control_mod, ONLY : iapp_tracvl,planet_type
    1212c
    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)

  • LMDZ5/trunk/libf/dyn3dpar/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/dyn3dpar/conf_gcm.F

    r1403 r1454  
    578578       offline = .FALSE.
    579579       CALL getin('offline',offline)
    580 
     580       IF (offline .AND. adjust) THEN
     581          WRITE(lunout,*)
     582     &         'WARNING : option offline does not work with adjust=y :'
     583          WRITE(lunout,*) 'the files defstoke.nc, fluxstoke.nc ',
     584     &         'and fluxstokev.nc will not be created'
     585          WRITE(lunout,*)
     586     &         'only the file phystoke.nc will still be created '
     587       END IF
     588       
    581589!Config  Key  = config_inca
    582590!Config  Desc = Choix de configuration de INCA
    … …  
    768776       offline = .FALSE.
    769777       CALL getin('offline',offline)
     778       IF (offline .AND. adjust) THEN
     779          WRITE(lunout,*)
     780     &         'WARNING : option offline does not work with adjust=y :'
     781          WRITE(lunout,*) 'the files defstoke.nc, fluxstoke.nc ',
     782     &         'and fluxstokev.nc will not be created'
     783          WRITE(lunout,*)
     784     &         'only the file phystoke.nc will still be created '
     785       END IF
    770786
    771787!Config  Key  = config_inca

  • LMDZ5/trunk/libf/dyn3dpar/ener.h

    r774 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/dyn3dpar/fluxstokenc_p.F

    r1403 r1454  
    5757
    5858c AC initialisations
    59 cym      pbarug(:,:)   = 0.
     59      pbarug(:,:)   = 0.
    6060cym      pbarvg(:,:,:) = 0.
    6161cym      wg(:,:)       = 0.

  • LMDZ5/trunk/libf/dyn3dpar/friction_p.F

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

    r1403 r1454  
    276276        endif
    277277
    278         if (planet_type.eq."earth") then
    279 #ifdef CPP_EARTH
     278!        if (planet_type.eq."earth") then
    280279! Load an Earth-format start file
    281280         CALL dynetat0("start.nc",vcov,ucov,
    282281     &              teta,q,masse,ps,phis, time_0)
    283 #else
    284         ! SW model also has Earth-format start files
    285         ! (but can be used without the CPP_EARTH directive)
    286           if (iflag_phys.eq.0) then
    287             CALL dynetat0("start.nc",vcov,ucov,
    288      &              teta,q,masse,ps,phis, time_0)
    289           endif
    290 #endif
    291         endif ! of if (planet_type.eq."earth")
     282!        endif ! of if (planet_type.eq."earth")
     283
    292284c       write(73,*) 'ucov',ucov
    293285c       write(74,*) 'vcov',vcov
    … …  
    494486#endif
    495487
    496       if (planet_type.eq."earth") then
     488!      if (planet_type.eq."earth") then
     489! Write an Earth-format restart file
    497490        CALL dynredem0_p("restart.nc", day_end, phis)
    498       endif
     491!      endif
    499492
    500493      ecripar = .TRUE.

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

    r1403 r1454  
    458458C MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
    459459
    460       PARAMETER (ISMo=300,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)'
    466463     
    467464      SUM=0.

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

    r1403 r1454  
    6161    CHARACTER(len=1), DIMENSION(3)  :: txts
    6262    CHARACTER(len=2), DIMENSION(9)  :: txtp
    63     CHARACTER(len=13)               :: str1,str2
     63    CHARACTER(len=23)               :: str1,str2
    6464 
    6565    INTEGER :: nqtrue  ! number of tracers read from tracer.def, without higer order of moment
    6666    INTEGER :: iq, new_iq, iiq, jq, ierr
    67     INTEGER, EXTERNAL :: lnblnk
    68  
     67
     68    character(len=*),parameter :: modname="infotrac_init"
    6969!-----------------------------------------------------------------------
    7070! Initialization :
    … …  
    100100       OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr)
    101101       IF(ierr.EQ.0) THEN
    102           WRITE(lunout,*) 'Open traceur.def : ok'
     102          WRITE(lunout,*) trim(modname),': Open traceur.def : ok'
    103103          READ(90,*) nqtrue
    104104       ELSE
    105           WRITE(lunout,*) 'Problem in opening traceur.def'
    106           WRITE(lunout,*) 'ATTENTION using defaut values'
    107           nqtrue=4 ! Defaut value
    108        END IF
    109        ! Attention! Only for planet_type=='earth'
    110        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
    111121    ELSE
    112122       ! nbtr has been read from INCA by init_cont_lmdz() in gcm.F
    … …  
    114124    END IF
    115125
    116     IF (nqtrue < 2) THEN
    117        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'
    118128       CALL abort_gcm('infotrac_init','Not enough tracers',1)
    119129    END IF
    … …  
    156166          ! Continue to read tracer.def
    157167          DO iq=1,nqtrue
    158              READ(90,999) hadv(iq),vadv(iq),tnom_0(iq)
     168             READ(90,*) hadv(iq),vadv(iq),tnom_0(iq)
    159169          END DO
    160170          CLOSE(90) 
    161        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
    162174          hadv(1) = 14
    163175          vadv(1) = 14
    … …  
    172184          vadv(4) = 10
    173185          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")
    174191       END IF
    175192       
    176        WRITE(lunout,*) 'Valeur de traceur.def :'
    177        WRITE(lunout,*) 'nombre de traceurs ',nqtrue
     193       WRITE(lunout,*) trim(modname),': Valeur de traceur.def :'
     194       WRITE(lunout,*) trim(modname),': nombre de traceurs ',nqtrue
    178195       DO iq=1,nqtrue
    179196          WRITE(lunout,*) hadv(iq),vadv(iq),tnom_0(iq)
    … …  
    217234          new_iq=new_iq+10 ! 9 tracers added
    218235       ELSE
    219           WRITE(lunout,*) 'This choice of advection schema is not available'
     236          WRITE(lunout,*) trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
    220237          CALL abort_gcm('infotrac_init','Bad choice of advection schema - 1',1)
    221238       END IF
    … …  
    227244       nqtot = new_iq
    228245
    229        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'
    230247       WRITE(lunout,*) 'makes it necessary to add tracers'
    231        WRITE(lunout,*) nqtrue,' is the number of true tracers'
    232        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'
    233250
    234251    ELSE
    … …  
    258275          iadv(new_iq)=11
    259276       ELSE
    260           WRITE(lunout,*)'This choice of advection schema is not available'
     277          WRITE(lunout,*)trim(modname),': This choice of advection schema is not available',iq,hadv(iq),vadv(iq)
     278
    261279          CALL abort_gcm('infotrac_init','Bad choice of advection schema - 2',1)
    262280       END IF
    … …  
    265283       tname(new_iq)= tnom_0(iq)
    266284       IF (iadv(new_iq)==0) THEN
    267           ttext(new_iq)=str1(1:lnblnk(str1))
     285          ttext(new_iq)=trim(str1)
    268286       ELSE
    269           ttext(new_iq)=str1(1:lnblnk(str1))//descrq(iadv(new_iq))
     287          ttext(new_iq)=trim(tnom_0(iq))//descrq(iadv(new_iq))
    270288       END IF
    271289
    … …  
    276294             new_iq=new_iq+1
    277295             iadv(new_iq)=-20
    278              ttext(new_iq)=str2(1:lnblnk(str2))//txts(jq)
    279              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)
    280298          END DO
    281299       ELSE IF (iadv(new_iq)==30) THEN
    … …  
    283301             new_iq=new_iq+1
    284302             iadv(new_iq)=-30
    285              ttext(new_iq)=str2(1:lnblnk(str2))//txtp(jq)
    286              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)
    287305          END DO
    288306       END IF
    … …  
    303321
    304322
    305     WRITE(lunout,*) 'Information stored in infotrac :'
    306     WRITE(lunout,*) 'iadv  niadv tname  ttext :'
     323    WRITE(lunout,*) trim(modname),': Information stored in infotrac :'
     324    WRITE(lunout,*) trim(modname),': iadv  niadv tname  ttext :'
    307325    DO iq=1,nqtot
    308        WRITE(lunout,*) iadv(iq),niadv(iq), tname(iq), ttext(iq)
     326       WRITE(lunout,*) iadv(iq),niadv(iq),&
     327       ' ',trim(tname(iq)),' ',trim(ttext(iq))
    309328    END DO
    310329
    … …  
    315334    DO iq=1,nqtot
    316335       IF (iadv(iq)/=10 .AND. iadv(iq)/=14 .AND. iadv(iq)/=0) THEN
    317           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'
    318337          CALL abort_gcm('infotrac_init','In this version only iadv=10 and iadv=14 is tested!',1)
    319338       ELSE IF (iadv(iq)==14 .AND. iq/=1) THEN
    320           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'
    321340          CALL abort_gcm('infotrac_init','In this version iadv=14 is only permitted for water vapour!',1)
    322341       END IF
    … …  
    329348    DEALLOCATE(tracnam)
    330349
    331 999 FORMAT (i2,1x,i2,1x,a15)
    332 
    333350  END SUBROUTINE infotrac_init
    334351

  • LMDZ5/trunk/libf/dyn3dpar/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/dyn3dpar/initfluxsto_p.F

    r1279 r1454  
    203203     .              llm, nivsigs, zvertiid)
    204204c pour le fichier def
    205       nivd(1) = 1
    206       call histvert(filedid, 'sig_s', 'Niveaux sigma',
    207      .  'sigma_level',
    208      .              1, nivd, dvertiid)
    209 
     205      if (mpi_rank==0) then
     206         nivd(1) = 1
     207         call histvert(filedid, 'sig_s', 'Niveaux sigma',
     208     .        'sigma_level',
     209     .        1, nivd, dvertiid)
     210      endif
    210211C
    211212C  Appels a histdef pour la definition des variables a sauvegarder
    … …  
    282283      call histend(fileid)
    283284      call histend(filevid)
    284       call histend(filedid)
     285      if (mpi_rank==0) call histend(filedid)
    285286      if (ok_sync) then
    286287        call histsync(fileid)
    287288        call histsync(filevid)
    288         call histsync(filedid)
     289        if (mpi_rank==0) call histsync(filedid)
    289290      endif
    290291       

  • LMDZ5/trunk/libf/dyn3dpar/integrd_p.F

    r1403 r1454  
    66     $     dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps0,masse,phis,finvmaold)
    77      USE parallel
    8       USE control_mod
     8      USE control_mod, only : planet_type
    99      IMPLICIT NONE
    1010
    … …  
    279279
    280280          CALL qminimum_p( q, nq, deltap )
    281          endif ! of if (planet_type.eq."earth")
    282281c
    283282c    .....  Calcul de la valeur moyenne, unique  aux poles pour  q .....
    … …  
    337336      ENDDO
    338337c$OMP END DO NOWAIT
     338
     339      endif ! of if (planet_type.eq."earth")
     340
    339341c
    340342c

  • LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F

    r1403 r1454  
    234234
    235235c$OMP MASTER
    236       dq=0.
     236      dq(:,:,:)=0.
    237237      CALL pression ( ip1jmp1, ap, bp, ps, p       )
    238238      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    … …  
    596596     .        flxw,pk, iapptrac)
    597597
    598        IF (offline) THEN
    599 Cmaf stokage du flux de masse pour traceurs OFF-LINE
    600 
    601 #ifdef CPP_IOIPSL
    602            CALL fluxstokenc_p(pbaru,pbarv,masse,teta,phi,phis,
    603      .   dtvr, itau)
    604 #endif
    605 
    606 
    607          ENDIF ! of IF (offline)
    608 c
     598C        Stokage du flux de masse pour traceurs OFF-LINE
     599         IF (offline .AND. .NOT. adjust) THEN
     600            CALL fluxstokenc_p(pbaru,pbarv,masse,teta,phi,phis,
     601     .           dtvr, itau)
     602         ENDIF
     603
    609604      ENDIF ! of IF( forward. OR . leapf )
    610605cc$OMP END PARALLEL
    … …  
    968963
    969964      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
     965!   Academic case : Simple friction and Newtonan relaxation
     966!   -------------------------------------------------------
    974967       ijb=ij_begin
    975968       ije=ij_end
    976969!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    977970       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
     971        teta(ijb:ije,l)=teta(ijb:ije,l)-dtvr*
     972     &         (teta(ijb:ije,l)-tetarappel(ijb:ije,l))*
     973     &                  (knewt_g+knewt_t(l)*clat4(ijb:ije))
     974       enddo ! of do l=1,llm
    981975!$OMP END DO
    982976
    … …  
    987981       call WaitRequest(Request_Physic)     
    988982c$OMP BARRIER
    989 !$OMP MASTER
    990        call friction_p(ucov,vcov,iphysiq*dtvr)
    991 !$OMP END MASTER
     983       call friction_p(ucov,vcov,dtvr)
    992984!$OMP BARRIER
     985
     986        ! Sponge layer (if any)
     987        IF (ok_strato) THEN
     988          ! set dufi,dvfi,... to zero
     989          ijb=ij_begin
     990          ije=ij_end
     991!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     992          do l=1,llm
     993            dufi(ijb:ije,l)=0
     994            dtetafi(ijb:ije,l)=0
     995            dqfi(ijb:ije,l,1:nqtot)=0
     996          enddo
     997!$OMP END DO
     998!$OMP SINGLE
     999          dpfi(ijb:ije)=0
     1000!$OMP END SINGLE
     1001          ijb=ij_begin
     1002          ije=ij_end
     1003          if (pole_sud) ije=ije-iip1
     1004!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     1005          do l=1,llm
     1006            dvfi(ijb:ije,l)=0
     1007          enddo
     1008!$OMP END DO
     1009
     1010          CALL top_bound_p(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
     1011          CALL addfi_p( dtvr, leapf, forward   ,
     1012     $                  ucov, vcov, teta , q   ,ps ,
     1013     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
     1014!$OMP BARRIER
     1015        ENDIF ! of IF (ok_strato)
    9931016      ENDIF ! of IF(iflag_phys.EQ.2)
    9941017
    … …  
    14651488c$OMP MASTER
    14661489
    1467               if (planet_type.eq."earth") then
     1490!              if (planet_type.eq."earth") then
    14681491! Write an Earth-format restart file
    14691492                CALL dynredem1_p("restart.nc",0.0,
    14701493     &                           vcov,ucov,teta,q,masse,ps)
    1471               endif ! of if (planet_type.eq."earth")
     1494!              endif ! of if (planet_type.eq."earth")
    14721495
    14731496!              CLOSE(99)
    … …  
    16581681
    16591682              IF(itau.EQ.itaufin) THEN
    1660                 if (planet_type.eq."earth") then
     1683!                if (planet_type.eq."earth") then
    16611684c$OMP MASTER
    16621685                   CALL dynredem1_p("restart.nc",0.0,
    16631686     .                               vcov,ucov,teta,q,masse,ps)
    16641687c$OMP END MASTER
    1665                 endif ! of if (planet_type.eq."earth")
     1688!                endif ! of if (planet_type.eq."earth")
    16661689              ENDIF ! of IF(itau.EQ.itaufin)
    16671690

  • LMDZ5/trunk/libf/dyn3dpar/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

  • LMDZ5/trunk/libf/phylmd/calcul_divers.h

    r1398 r1454  
    22c $Header$
    33c
    4 c
    5 c initialisations diverses au "debut" du mois
    6 c
    7       IF(MOD(itap,INT(ecrit_mth/dtime)).EQ.1) THEN
     4
     5c     Initialisations diverses au "debut" du mois
     6      IF(debut) THEN
     7         nday_rain(:)=0.
     8
     9c        surface terre
     10         paire_ter(:)=0.
    811         DO i=1, klon
    9           nday_rain(i)=0.
    10          ENDDO !i
    11 c
    12 c surface terre
    13        DO i=1, klon
    14          IF(pctsrf(i,is_ter).GT.0.) THEN
    15             paire_ter(i)=airephy(i)*pctsrf(i,is_ter)
    16          ENDIF
    17        ENDDO
    18 c
    19       ENDIF !MOD(itap,INT(ecrit_mth)).EQ.1
    20 c
    21       IF(MOD(itap,INT(ecrit_day/dtime)).EQ.0) THEN
    22 c
    23 cIM calcul total_rain, nday_rain
    24 c
    25        DO i = 1, klon
    26         total_rain(i)=rain_fall(i)+snow_fall(i) 
    27         IF(total_rain(i).GT.0.) nday_rain(i)=nday_rain(i)+1.
    28        ENDDO
    29       ENDIF !itap.EQ.ecrit_mth
     12            IF(pctsrf(i,is_ter).GT.0.) THEN
     13               paire_ter(i)=airephy(i)*pctsrf(i,is_ter)
     14            ENDIF
     15         ENDDO
     16      ENDIF
     17
     18cIM   Calcul une fois par jour : total_rain, nday_rain
     19      IF(MOD(itap,INT(un_jour/dtime)).EQ.0) THEN
     20         DO i = 1, klon
     21            total_rain(i)=rain_fall(i)+snow_fall(i) 
     22            IF(total_rain(i).GT.0.) nday_rain(i)=nday_rain(i)+1.
     23         ENDDO
     24      ENDIF

  • LMDZ5/trunk/libf/phylmd/carbon_cycle_mod.F90

    r1279 r1454  
    11MODULE carbon_cycle_mod
    2 
     2! Controle module for the carbon CO2 tracers :
     3!   - Identification
     4!   - Get concentrations comming from coupled model or read from file to tracers
     5!   - Calculate new RCO2 for radiation scheme
     6!   - Calculate new carbon flux for sending to coupled models (PISCES and ORCHIDEE)
     7!
    38! Author : Josefine GHATTAS, Patricia CADULE
    49
    … …  
    1318  LOGICAL, PUBLIC :: carbon_cycle_cpl       ! Coupling of CO2 fluxes between LMDZ/ORCHIDEE and LMDZ/OCEAN(PISCES)
    1419!$OMP THREADPRIVATE(carbon_cycle_cpl)
     20
    1521  LOGICAL :: carbon_cycle_emis_comp=.FALSE. ! Calculation of emission compatible
     22!$OMP THREADPRIVATE(carbon_cycle_emis_comp)
     23
     24  LOGICAL :: RCO2_inter  ! RCO2 interactive : if true calculate new value RCO2 for the radiation scheme
     25!$OMP THREADPRIVATE(RCO2_inter)
    1626
    1727! Scalare values when no transport, from physiq.def
    … …  
    2131!$OMP THREADPRIVATE(emis_land_s)
    2232
    23   INTEGER :: ntr_co2                ! Number of tracers concerning the carbon cycle
    24   INTEGER :: id_fco2_tot            ! Tracer index
    25   INTEGER :: id_fco2_ocn            !  - " -
    26   INTEGER :: id_fco2_land           !  - " -
    27   INTEGER :: id_fco2_land_use       !  - " -
    28   INTEGER :: id_fco2_fos_fuel       !  - " -
    29 !$OMP THREADPRIVATE(ntr_co2, id_fco2_tot, id_fco2_ocn, id_fco2_land, id_fco2_land_use, id_fco2_fos_fuel) 
    30 
    31   REAL, DIMENSION(:), ALLOCATABLE :: fos_fuel        ! CO2 fossil fuel emission from file [gC/m2/d]
    32 !$OMP THREADPRIVATE(fos_fuel)
    33   REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocn_day ! flux CO2 from ocean for 1 day (cumulated) [gC/m2/d]
     33  REAL :: airetot     ! Total area of the earth surface
     34!$OMP THREADPRIVATE(airetot)
     35
     36  INTEGER :: ntr_co2  ! Number of tracers concerning the carbon cycle
     37!$OMP THREADPRIVATE(ntr_co2)
     38
     39! fco2_ocn_day : flux CO2 from ocean for 1 day (cumulated) [gC/m2/d]. Allocation and initalization done in cpl_mod
     40  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_ocn_day
    3441!$OMP THREADPRIVATE(fco2_ocn_day)
     42
    3543  REAL, DIMENSION(:), ALLOCATABLE :: fco2_land_day   ! flux CO2 from land for 1 day (cumulated)  [gC/m2/d]
    3644!$OMP THREADPRIVATE(fco2_land_day)
    … …  
    3846!$OMP THREADPRIVATE(fco2_lu_day)
    3947
    40 ! Following 2 fields will be initialized in surf_land_orchidee at each time step
     48  REAL, DIMENSION(:,:), ALLOCATABLE :: dtr_add       ! Tracer concentration to be injected
     49!$OMP THREADPRIVATE(dtr_add)
     50
     51! Following 2 fields will be allocated and initialized in surf_land_orchidee
    4152  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: fco2_land_inst  ! flux CO2 from land at one time step
    4253!$OMP THREADPRIVATE(fco2_land_inst)
    … …  
    4556
    4657! Calculated co2 field to be send to the ocean via the coupler and to ORCHIDEE
    47   REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: co2_send
     58  REAL, DIMENSION(:), ALLOCATABLE, PUBLIC :: co2_send ! Field allocated in phyetat0
    4859!$OMP THREADPRIVATE(co2_send)
     60
     61
     62  TYPE, PUBLIC ::   co2_trac_type
     63     CHARACTER(len = 8) :: name       ! Tracer name in tracer.def
     64     INTEGER            :: id         ! Index in total tracer list, tr_seri
     65     CHARACTER(len=30)  :: file       ! File name
     66     LOGICAL            :: cpl        ! True if this tracers is coupled from ORCHIDEE or PISCES.
     67                                      ! False if read from file.
     68     INTEGER            :: updatefreq ! Frequence to inject in second
     69     INTEGER            :: readstep   ! Actual time step to read in file
     70     LOGICAL            :: updatenow  ! True if this tracer should be updated this time step
     71  END TYPE co2_trac_type
     72  INTEGER,PARAMETER :: maxco2trac=5  ! Maximum number of different CO2 fluxes
     73  TYPE(co2_trac_type), DIMENSION(maxco2trac) :: co2trac
    4974
    5075CONTAINS
    5176 
    52   SUBROUTINE carbon_cycle_init(tr_seri, aerosol, radio)
     77  SUBROUTINE carbon_cycle_init(tr_seri, pdtphys, aerosol, radio)
     78! This subroutine is called from traclmdz_init, only at first timestep.
     79! - Read controle parameters from .def input file
     80! - Search for carbon tracers and set default values
     81! - Allocate variables
     82! - Test for compatibility
     83
    5384    USE dimphy
     85    USE comgeomphy
     86    USE mod_phys_lmdz_transfert_para
    5487    USE infotrac
    5588    USE IOIPSL
    5689    USE surface_data, ONLY : ok_veget, type_ocean
     90    USE phys_cal_mod, ONLY : mth_len
    5791
    5892    IMPLICIT NONE
    … …  
    6296! Input argument
    6397    REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr_seri ! Concentration Traceur [U/KgA] 
     98    REAL,INTENT(IN)                           :: pdtphys ! length of time step in physiq (sec)
    6499
    65100! InOutput arguments
    … …  
    68103
    69104! Local variables
    70     INTEGER               :: ierr, it, iiq
    71     REAL, DIMENSION(klon) :: tr_seri_sum
    72 
    73 
    74 ! 0) Test for compatibility
    75     IF (carbon_cycle_cpl .AND. type_ocean/='couple') &
    76          CALL abort_gcm('carbon_cycle_init', 'Coupling with ocean model is needed for carbon_cycle_cpl',1)
    77     IF (carbon_cycle_cpl .AND..NOT. ok_veget) &
    78          CALL abort_gcm('carbon_cycle_init', 'Coupling with surface land model ORCHDIEE is needed for carbon_cycle_cpl',1)
    79 
    80 
    81 ! 1) Check if transport of one tracer flux CO2 or 4 separated tracers
    82     IF (carbon_cycle_tr) THEN
    83        id_fco2_tot=0
    84        id_fco2_ocn=0
    85        id_fco2_land=0
    86        id_fco2_land_use=0
    87        id_fco2_fos_fuel=0
    88        
    89        ! Search in tracer list
    90        DO it=1,nbtr
    91           iiq=niadv(it+2)
    92           IF (tname(iiq) == "fCO2" ) THEN
    93              id_fco2_tot=it
    94           ELSE IF (tname(iiq) == "fCO2_ocn" ) THEN
    95              id_fco2_ocn=it
    96           ELSE IF (tname(iiq) == "fCO2_land" ) THEN
    97              id_fco2_land=it
    98           ELSE IF (tname(iiq) == "fCO2_land_use" ) THEN
    99              id_fco2_land_use=it
    100           ELSE IF (tname(iiq) == "fCO2_fos_fuel" ) THEN
    101              id_fco2_fos_fuel=it
    102           END IF
    103        END DO
    104 
    105        ! Count tracers found
    106        IF (id_fco2_tot /= 0 .AND. &
    107             id_fco2_ocn==0 .AND. id_fco2_land==0 .AND. id_fco2_land_use==0 .AND. id_fco2_fos_fuel==0) THEN
    108          
    109           ! transport  1 tracer flux CO2
    110           ntr_co2 = 1
    111          
    112        ELSE IF (id_fco2_tot==0 .AND. &
    113             id_fco2_ocn /=0 .AND. id_fco2_land/=0 .AND. id_fco2_land_use/=0 .AND. id_fco2_fos_fuel/=0) THEN
    114          
    115           ! transport 4 tracers seperatively
    116           ntr_co2 = 4
    117          
    118        ELSE
    119           CALL abort_gcm('carbon_cycle_init', 'error in coherence between traceur.def and gcm.def',1)
    120        END IF
    121 
    122        ! Definition of control varaiables for the tracers
    123        DO it=1,nbtr
    124           IF (it==id_fco2_tot .OR. it==id_fco2_ocn .OR. it==id_fco2_land .OR. &
    125                it==id_fco2_land_use .OR. it==id_fco2_fos_fuel) THEN
    126              aerosol(it) = .FALSE.
    127              radio(it)   = .FALSE.
    128           END IF
    129        END DO
    130 
    131     ELSE
    132        ! No transport of CO2
    133        ntr_co2 = 0
    134     END IF ! carbon_cycle_tr
    135 
    136 
    137 ! 2) Allocate variable for CO2 fossil fuel emission
    138     IF (carbon_cycle_tr) THEN
    139        ! Allocate 2D variable
    140        ALLOCATE(fos_fuel(klon), stat=ierr)
    141        IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 1',1)
    142     ELSE
    143        ! No transport : read value from .def
     105    INTEGER               :: ierr, it, iiq, itc
     106    INTEGER               :: teststop
     107
     108
     109
     110! 1) Read controle parameters from .def input file
     111! ------------------------------------------------
     112    ! Read fosil fuel value if no transport
     113    IF (.NOT. carbon_cycle_tr) THEN
    144114       fos_fuel_s = 0.
    145115       CALL getin ('carbon_cycle_fos_fuel',fos_fuel_s)
    … …  
    148118
    149119
    150 ! 3) Allocate and initialize fluxes
    151     IF (carbon_cycle_cpl) THEN
    152        IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 2',1)
    153        ALLOCATE(fco2_land_day(klon), stat=ierr)
    154        IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 3',1)
    155        ALLOCATE(fco2_lu_day(klon), stat=ierr)
    156        IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 4',1)
    157 
    158        fco2_land_day(:) = 0.  ! JG : Doit prend valeur de restart
    159        fco2_lu_day(:)   = 0.  ! JG : Doit prend valeur de restart
    160 
    161        ! fco2_ocn_day is allocated in cpl_mod
    162        ! fco2_land_inst and fco2_lu_inst are allocated in surf_land_orchidee
    163        
    164        ALLOCATE(co2_send(klon), stat=ierr)
    165        IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 7',1)
    166        
    167        ! Calculate using restart tracer values
    168        IF (carbon_cycle_tr) THEN
    169           IF (ntr_co2==1) THEN
    170              co2_send(:) = tr_seri(:,1,id_fco2_tot) + co2_ppm0
    171           ELSE ! ntr_co2==4
    172              ! Calculate the delta CO2 flux
    173              tr_seri_sum(:) = tr_seri(:,1,id_fco2_fos_fuel) + tr_seri(:,1,id_fco2_land_use) + &
    174                   tr_seri(:,1,id_fco2_land) + tr_seri(:,1,id_fco2_ocn)
    175              co2_send(:) = tr_seri_sum(:) + co2_ppm0
    176           END IF
    177        ELSE
    178           ! Send a scalare value in 2D variable to ocean and land model (PISCES and ORCHIDEE)
    179           co2_send(:) = co2_ppm
    180        END IF
    181 
    182 
    183     ELSE
    184        IF (carbon_cycle_tr) THEN
    185           ! No coupling of CO2 fields :
    186           ! corresponding fields will instead be read from files
    187           ALLOCATE(fco2_ocn_day(klon), stat=ierr)
    188           IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 8',1)
    189           ALLOCATE(fco2_land_day(klon), stat=ierr)
    190           IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 9',1)
    191           ALLOCATE(fco2_lu_day(klon), stat=ierr)
    192           IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 10',1)       
    193        END IF
    194     END IF
    195 
    196 ! 4) Read parmeter for calculation of emission compatible
     120    ! Read parmeter for calculation compatible emission
    197121    IF (.NOT. carbon_cycle_tr) THEN
    198122       carbon_cycle_emis_comp=.FALSE.
    199123       CALL getin('carbon_cycle_emis_comp',carbon_cycle_emis_comp)
    200124       WRITE(lunout,*) 'carbon_cycle_emis_comp = ',carbon_cycle_emis_comp
    201     END IF
     125       IF (carbon_cycle_emis_comp) THEN
     126          CALL abort_gcm('carbon_cycle_init', 'carbon_cycle_emis_comp option not yet implemented!!',1)
     127       END IF
     128    END IF
     129
     130    ! Read parameter for interactive calculation of the CO2 value for the radiation scheme
     131    RCO2_inter=.FALSE.
     132    CALL getin('RCO2_inter',RCO2_inter)
     133    WRITE(lunout,*) 'RCO2_inter = ', RCO2_inter
     134    IF (RCO2_inter) THEN
     135       WRITE(lunout,*) 'RCO2 will be recalculated once a day'
     136       WRITE(lunout,*) 'RCO2 initial = ', RCO2
     137    END IF
     138
     139
     140! 2) Search for carbon tracers and set default values
     141! ---------------------------------------------------
     142    itc=0
     143    DO it=1,nbtr
     144       iiq=niadv(it+2)
     145       
     146       SELECT CASE(tname(iiq))
     147       CASE("fCO2_ocn")
     148          itc = itc + 1
     149          co2trac(itc)%name='fCO2_ocn'
     150          co2trac(itc)%id=it
     151          co2trac(itc)%file='fl_co2_ocean.nc'
     152          IF (carbon_cycle_cpl .AND. type_ocean=='couple') THEN
     153             co2trac(itc)%cpl=.TRUE.
     154             co2trac(itc)%updatefreq = 86400 ! Once a day as the coupling with OASIS/PISCES
     155          ELSE
     156             co2trac(itc)%cpl=.FALSE.
     157             co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
     158          END IF
     159       CASE("fCO2_land")
     160          itc = itc + 1
     161          co2trac(itc)%name='fCO2_land'
     162          co2trac(itc)%id=it
     163          co2trac(itc)%file='fl_co2_land.nc'
     164          IF (carbon_cycle_cpl .AND. ok_veget) THEN
     165             co2trac(itc)%cpl=.TRUE.
     166             co2trac(itc)%updatefreq = INT(pdtphys) ! Each timestep as the coupling with ORCHIDEE
     167          ELSE
     168             co2trac(itc)%cpl=.FALSE.
     169!             co2trac(itc)%updatefreq = 10800   ! 10800sec = 3H
     170             co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
     171          END IF
     172       CASE("fCO2_land_use")
     173          itc = itc + 1
     174          co2trac(itc)%name='fCO2_land_use'
     175          co2trac(itc)%id=it
     176          co2trac(itc)%file='fl_co2_land_use.nc'
     177          IF (carbon_cycle_cpl .AND. ok_veget) THEN
     178             co2trac(it)%cpl=.TRUE.
     179             co2trac(itc)%updatefreq = INT(pdtphys) ! Each timestep as the coupling with ORCHIDEE
     180          ELSE
     181             co2trac(itc)%cpl=.FALSE.
     182             co2trac(itc)%updatefreq = 10800   ! 10800sec = 3H
     183          END IF
     184       CASE("fCO2_fos_fuel")
     185          itc = itc + 1
     186          co2trac(itc)%name='fCO2_fos_fuel'
     187          co2trac(itc)%id=it
     188          co2trac(itc)%file='fossil_fuel.nc'
     189          co2trac(itc)%cpl=.FALSE.       ! This tracer always read from file
     190!         co2trac(itc)%updatefreq = 86400  ! 86400sec = 24H Cadule case
     191          co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
     192       CASE("fCO2_bbg")
     193          itc = itc + 1
     194          co2trac(itc)%name='fCO2_bbg'
     195          co2trac(itc)%id=it
     196          co2trac(itc)%file='fl_co2_bbg.nc'
     197          co2trac(itc)%cpl=.FALSE.       ! This tracer always read from file
     198          co2trac(itc)%updatefreq = 86400*mth_len ! Once a month
     199       CASE("fCO2")
     200          ! fCO2 : One tracer transporting the total CO2 flux
     201          itc = itc + 1
     202          co2trac(itc)%name='fCO2'
     203          co2trac(itc)%id=it
     204          co2trac(itc)%file='fl_co2.nc'
     205          IF (carbon_cycle_cpl) THEN
     206             co2trac(itc)%cpl=.TRUE.
     207          ELSE
     208             co2trac(itc)%cpl=.FALSE.
     209          END IF
     210          co2trac(itc)%updatefreq = 86400
     211          ! DOES THIS WORK ???? Problematic due to implementation of the coupled fluxes...
     212          CALL abort_gcm('carbon_cycle_init','transport of total CO2 has to be implemented and tested',1)
     213       END SELECT
     214    END DO
     215
     216    ! Total number of carbon CO2 tracers
     217    ntr_co2 = itc
     218   
     219    ! Definition of control varaiables for the tracers
     220    DO it=1,ntr_co2
     221       aerosol(co2trac(it)%id) = .FALSE.
     222       radio(co2trac(it)%id)   = .FALSE.
     223    END DO
     224   
     225    ! Vector indicating which timestep to read for each tracer
     226    ! Always start read in the beginning of the file
     227    co2trac(:)%readstep = 0
     228   
     229
     230! 3) Allocate variables
     231! ---------------------
     232    ! Allocate vector for storing fluxes to inject
     233    ALLOCATE(dtr_add(klon,maxco2trac), stat=ierr)
     234    IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 11',1)       
     235   
     236    ! Allocate variables for cumulating fluxes from ORCHIDEE
     237    IF (RCO2_inter) THEN
     238       IF (.NOT. carbon_cycle_tr .AND. carbon_cycle_cpl) THEN
     239          ALLOCATE(fco2_land_day(klon), stat=ierr)
     240          IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 2',1)
     241          fco2_land_day(1:klon) = 0.
     242         
     243          ALLOCATE(fco2_lu_day(klon), stat=ierr)
     244          IF (ierr /= 0) CALL abort_gcm('carbon_cycle_init', 'pb in allocation 3',1)
     245          fco2_lu_day(1:klon)   = 0.
     246       END IF
     247    END IF
     248
     249
     250! 4) Test for compatibility
     251! -------------------------
     252!    IF (carbon_cycle_cpl .AND. type_ocean/='couple') THEN
     253!       WRITE(lunout,*) 'Coupling with ocean model is needed for carbon_cycle_cpl'
     254!       CALL abort_gcm('carbon_cycle_init', 'coupled ocean is needed for carbon_cycle_cpl',1)
     255!    END IF
     256!
     257!    IF (carbon_cycle_cpl .AND..NOT. ok_veget) THEN
     258!       WRITE(lunout,*) 'Coupling with surface land model ORCHDIEE is needed for carbon_cycle_cpl'
     259!       CALL abort_gcm('carbon_cycle_init', 'ok_veget is needed for carbon_cycle_cpl',1)
     260!    END IF
     261
     262    ! Compiler test : following should never happen
     263    teststop=0
     264    DO it=1,teststop
     265       CALL abort_gcm('carbon_cycle_init', 'Entering loop from 1 to 0',1)
     266    END DO
     267
     268    IF (ntr_co2==0) THEN
     269       ! No carbon tracers found in tracer.def. It is not possible to do carbon cycle
     270       WRITE(lunout,*) 'No carbon tracers found in tracer.def. Not ok with carbon_cycle_tr and/or carbon_cycle_cp'
     271       CALL abort_gcm('carbon_cycle_init', 'No carbon tracers found in tracer.def',1)
     272    END IF
     273   
     274! 5) Calculate total area of the earth surface
     275! --------------------------------------------
     276    CALL reduce_sum(SUM(airephy),airetot)
     277    CALL bcast(airetot)
    202278
    203279  END SUBROUTINE carbon_cycle_init
    204280
     281
     282  SUBROUTINE carbon_cycle(nstep, pdtphys, pctsrf, tr_seri, source)
     283! Subroutine for injection of co2 in the tracers
    205284!
    206 !
    207 !
    208 
    209   SUBROUTINE carbon_cycle(nstep, pdtphys, pctsrf, tr_seri)
    210    
     285! - Find out if it is time to update
     286! - Get tracer from coupled model or from file
     287! - Calculate new RCO2 value for the radiation scheme
     288! - Calculate CO2 flux to send to ocean and land models (PISCES and ORCHIDEE)
     289
    211290    USE infotrac
    212291    USE dimphy
    213     USE mod_phys_lmdz_transfert_para, ONLY : reduce_sum
     292    USE mod_phys_lmdz_transfert_para
    214293    USE phys_cal_mod, ONLY : mth_cur, mth_len
    215294    USE phys_cal_mod, ONLY : day_cur
    … …  
    220299    INCLUDE "clesphys.h"
    221300    INCLUDE "indicesol.h"
     301    INCLUDE "iniprint.h"
     302    INCLUDE "YOMCST.h"
    222303
    223304! In/Output arguments
    224305    INTEGER,INTENT(IN) :: nstep      ! time step in physiq
    225306    REAL,INTENT(IN)    :: pdtphys    ! length of time step in physiq (sec)
    226     REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
    227     REAL, DIMENSION(klon,klev,nbtr), INTENT(INOUT)  :: tr_seri
     307    REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf            ! Surface fraction
     308    REAL, DIMENSION(klon,klev,nbtr), INTENT(INOUT)  :: tr_seri ! All tracers
     309    REAL, DIMENSION(klon,nbtr), INTENT(INOUT)       :: source  ! Source for all tracers
    228310
    229311! Local variables
     312    INTEGER :: it
    230313    LOGICAL :: newmonth ! indicates if a new month just started
    231314    LOGICAL :: newday   ! indicates if a new day just started
    … …  
    233316
    234317    REAL, PARAMETER :: fact=1.E-15/2.12  ! transformation factor from gC/m2/day => ppm/m2/day
    235     REAL, DIMENSION(klon) :: fco2_tmp, tr_seri_sum
     318    REAL, DIMENSION(klon) :: fco2_tmp
    236319    REAL :: sumtmp
    237     REAL :: airetot     ! Total area the earth
    238320    REAL :: delta_co2_ppm
    239321   
    240 ! -) Calculate logicals indicating if it is a new month, new day or the last time step in a day (end day)
     322
     323! 1) Calculate logicals indicating if it is a new month, new day or the last time step in a day (end day)
     324! -------------------------------------------------------------------------------------------------------
    241325
    242326    newday = .FALSE.; endday = .FALSE.; newmonth = .FALSE.
    … …  
    245329    IF (MOD(nstep,INT(86400./pdtphys))==0) endday=.TRUE.
    246330    IF (newday .AND. day_cur==1) newmonth=.TRUE.
    247    
    248 ! -) Read new maps if new month started
    249     IF (newmonth .AND. carbon_cycle_tr) THEN
    250        CALL read_map2D('fossil_fuel.nc','fos_fuel', mth_cur, .FALSE., fos_fuel)
    251        
    252        ! division by month lenght to get dayly value
    253        fos_fuel(:) = fos_fuel(:)/mth_len
    254        
    255        IF (.NOT. carbon_cycle_cpl) THEN
    256           ! Get dayly values from monthly fluxes
    257           CALL read_map2D('fl_co2_ocean.nc','CO2_OCN',mth_cur,.FALSE.,fco2_ocn_day)
    258           CALL read_map2D('fl_co2_land.nc','CO2_LAND', mth_cur,.FALSE.,fco2_land_day)
    259           CALL read_map2D('fl_co2_land_use.nc','CO2_LAND_USE',mth_cur,.FALSE.,fco2_lu_day)
    260        END IF
    261     END IF
    262    
    263 
    264 
    265 ! -) Update tracers at beginning of a new day. Beginning of a new day correspond to a new coupling period in cpl_mod.
    266     IF (newday) THEN
     331
     332! 2)  For each carbon tracer find out if it is time to inject (update)
     333! --------------------------------------------------------------------
     334    DO it = 1, ntr_co2
     335       IF ( MOD(nstep,INT(co2trac(it)%updatefreq/pdtphys)) == 1 ) THEN
     336          co2trac(it)%updatenow = .TRUE.
     337       ELSE
     338          co2trac(it)%updatenow = .FALSE.
     339       END IF
     340    END DO
     341
     342! 3) Get tracer update
     343! --------------------------------------
     344    DO it = 1, ntr_co2
     345       IF ( co2trac(it)%updatenow ) THEN
     346          IF ( co2trac(it)%cpl ) THEN
     347             ! Get tracer from coupled model
     348             SELECT CASE(co2trac(it)%name)
     349             CASE('fCO2_land')     ! from ORCHIDEE
     350                dtr_add(:,it) = fco2_land_inst(:)*pctsrf(:,is_ter)*fact ! [ppm/m2/day]
     351             CASE('fCO2_land_use') ! from ORCHIDEE
     352                dtr_add(:,it) = fco2_lu_inst(:)  *pctsrf(:,is_ter)*fact ! [ppm/m2/day]
     353             CASE('fCO2_ocn')      ! from PISCES
     354                dtr_add(:,it) = fco2_ocn_day(:)  *pctsrf(:,is_oce)*fact ! [ppm/m2/day]
     355             CASE DEFAULT
     356                WRITE(lunout,*) 'Error with tracer ',co2trac(it)%name
     357                CALL abort_gcm('carbon_cycle', 'No coupling implemented for this tracer',1)
     358             END SELECT
     359          ELSE
     360             ! Read tracer from file
     361             co2trac(it)%readstep = co2trac(it)%readstep + 1 ! increment time step in file
     362! Patricia   CALL read_map2D(co2trac(it)%file,'fco2',co2trac(it)%readstep,.FALSE.,dtr_add(:,it))
     363             CALL read_map2D(co2trac(it)%file,'fco2',co2trac(it)%readstep,.TRUE.,dtr_add(:,it))
     364
     365             ! Converte from kgC/m2/h to kgC/m2/s
     366             dtr_add(:,it) = dtr_add(:,it)/3600
     367             ! Add individual treatment of values read from file
     368             SELECT CASE(co2trac(it)%name)
     369             CASE('fCO2_land')
     370                dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_ter)
     371             CASE('fCO2_land_use')
     372                dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_ter)
     373             CASE('fCO2_ocn')
     374                dtr_add(:,it) = dtr_add(:,it) *pctsrf(:,is_oce)
     375! Patricia :
     376!             CASE('fCO2_fos_fuel')
     377!                dtr_add(:,it) = dtr_add(:,it)/mth_len
     378!                co2trac(it)%readstep = 0 ! Always read same value for fossil fuel(Cadule case)
     379             END SELECT
     380          END IF
     381       END IF
     382    END DO
     383
     384! 4) Update co2 tracers :
     385!    Loop over all carbon tracers and add source
     386! ------------------------------------------------------------------
     387    IF (carbon_cycle_tr) THEN
     388       DO it = 1, ntr_co2
     389          IF (.FALSE.) THEN
     390             tr_seri(1:klon,1,co2trac(it)%id) = tr_seri(1:klon,1,co2trac(it)%id) + dtr_add(1:klon,it)
     391             source(1:klon,co2trac(it)%id) = 0.
     392          ELSE
     393             source(1:klon,co2trac(it)%id) = dtr_add(1:klon,it)
     394          END IF
     395       END DO
     396    END IF
     397
     398
     399! 5) Calculations for new CO2 value for the radiation scheme(instead of reading value from .def)
     400! ----------------------------------------------------------------------------------------------
     401    IF (RCO2_inter) THEN
     402       ! Cumulate fluxes from ORCHIDEE at each timestep
     403       IF (.NOT. carbon_cycle_tr .AND. carbon_cycle_cpl) THEN
     404          IF (newday) THEN ! Reset cumulative variables once a day
     405             fco2_land_day(1:klon) = 0.
     406             fco2_lu_day(1:klon)   = 0.
     407          END IF
     408          fco2_land_day(1:klon) = fco2_land_day(1:klon) + fco2_land_inst(1:klon) ![gC/m2/day]
     409          fco2_lu_day(1:klon)   = fco2_lu_day(1:klon)   + fco2_lu_inst(1:klon)   ![gC/m2/day]
     410       END IF
     411
     412       ! At the end of a new day, calculate a mean scalare value of CO2
     413       ! JG : Ici on utilise uniquement le traceur du premier couche du modele. Est-ce que c'est correcte ?
     414       IF (endday) THEN
     415
     416          IF (carbon_cycle_tr) THEN
     417             ! Sum all co2 tracers to get the total delta CO2 flux
     418             fco2_tmp(:) = 0.
     419             DO it = 1, ntr_co2
     420                fco2_tmp(1:klon) = fco2_tmp(1:klon) + tr_seri(1:klon,1,co2trac(it)%id)
     421             END DO
     422             
     423          ELSE IF (carbon_cycle_cpl) THEN ! no carbon_cycle_tr
     424             ! Sum co2 fluxes comming from coupled models and parameter for fossil fuel
     425             fco2_tmp(1:klon) = fos_fuel_s + ((fco2_lu_day(1:klon) + fco2_land_day(1:klon))*pctsrf(1:klon,is_ter) &
     426                  + fco2_ocn_day(:)*pctsrf(:,is_oce)) * fact
     427          END IF
     428
     429          ! Calculate a global mean value of delta CO2 flux
     430          fco2_tmp(1:klon) = fco2_tmp(1:klon) * airephy(1:klon)
     431          CALL reduce_sum(SUM(fco2_tmp),sumtmp)
     432          CALL bcast(sumtmp)
     433          delta_co2_ppm = sumtmp/airetot
     434         
     435          ! Add initial value for co2_ppm and delta value
     436          co2_ppm = co2_ppm0 + delta_co2_ppm
     437         
     438          ! Transformation of atmospheric CO2 concentration for the radiation code
     439          RCO2 = co2_ppm * 1.0e-06  * 44.011/28.97
     440         
     441          WRITE(lunout,*) 'RCO2 is now updated! RCO2 = ', RCO2
     442       END IF ! endday
     443
     444    END IF ! RCO2_inter
     445
     446
     447! 6) Calculate CO2 flux to send to ocean and land models : PISCES and ORCHIDEE         
     448! ----------------------------------------------------------------------------
     449    IF (carbon_cycle_cpl) THEN
    267450
    268451       IF (carbon_cycle_tr) THEN
    269 
    270           ! Update tracers
    271           IF (ntr_co2 == 1) THEN
    272              ! Calculate the new flux CO2
    273              tr_seri(:,1,id_fco2_tot) = tr_seri(:,1,id_fco2_tot) + &
    274                   (fos_fuel(:) + &
    275                   fco2_lu_day(:)  * pctsrf(:,is_ter) + &
    276                   fco2_land_day(:)* pctsrf(:,is_ter) + &
    277                   fco2_ocn_day(:) * pctsrf(:,is_oce)) * fact
    278 
    279           ELSE ! ntr_co2 == 4
    280              tr_seri(:,1,id_fco2_fos_fuel)  = tr_seri(:,1,id_fco2_fos_fuel) + fos_fuel(:) * fact ! [ppm/m2/day]
    281 
    282              tr_seri(:,1,id_fco2_land_use)  = tr_seri(:,1,id_fco2_land_use) + &
    283                   fco2_lu_day(:)  *pctsrf(:,is_ter)*fact ! [ppm/m2/day]
    284 
    285              tr_seri(:,1,id_fco2_land)      = tr_seri(:,1,id_fco2_land)     + &
    286                   fco2_land_day(:)*pctsrf(:,is_ter)*fact ! [ppm/m2/day]
    287 
    288              tr_seri(:,1,id_fco2_ocn)       = tr_seri(:,1,id_fco2_ocn)      + &
    289                   fco2_ocn_day(:) *pctsrf(:,is_oce)*fact ! [ppm/m2/day]
    290 
    291           END IF
    292 
    293        ELSE ! no transport
    294           IF (carbon_cycle_cpl) THEN
    295              IF (carbon_cycle_emis_comp) THEN
    296                 ! Calcul emission compatible a partir des champs 2D et co2_ppm
    297                 !!! TO DO!!
    298                 CALL abort_gcm('carbon_cycle', ' Option carbon_cycle_emis_comp not yet implemented',1)
    299              END IF
    300           END IF
    301 
    302        END IF ! carbon_cycle_tr
    303    
    304        ! Reset cumluative variables
    305        IF (carbon_cycle_cpl) THEN
    306           fco2_land_day(:) = 0.
    307           fco2_lu_day(:)   = 0.
    308        END IF
    309    
    310     END IF ! newday
    311        
    312 
    313 
    314 ! -) Cumulate fluxes from ORCHIDEE at each timestep
    315     IF (carbon_cycle_cpl) THEN
    316        fco2_land_day(:) = fco2_land_day(:) + fco2_land_inst(:)
    317        fco2_lu_day(:)   = fco2_lu_day(:)   + fco2_lu_inst(:)
    318     END IF
    319 
    320 
    321 
    322 ! -)  At the end of a new day, calculate a mean scalare value of CO2 to be used by
    323 !     the radiation scheme (instead of reading value from .def)
    324 
    325 ! JG : Ici on utilise uniquement le traceur du premier couche du modele. Est-ce que c'est correcte ?
    326     IF (endday) THEN
    327        ! Calculte total area of the earth surface
    328        CALL reduce_sum(SUM(airephy),airetot)
    329        
    330 
    331        IF (carbon_cycle_tr) THEN
    332           IF (ntr_co2 == 1) THEN
    333          
    334              ! Calculate mean value of tracer CO2 to get an scalare value to be used in the
    335              ! radiation scheme (instead of reading value from .def)
    336              ! Mean value weighted with the grid cell area
    337              
    338              ! Calculate mean value
    339              fco2_tmp(:) = tr_seri(:,1,id_fco2_tot) * airephy(:)
    340              CALL reduce_sum(SUM(fco2_tmp),sumtmp)
    341              co2_ppm = sumtmp/airetot + co2_ppm0
    342              
    343           ELSE ! ntr_co2 == 4
    344              
    345              ! Calculate the delta CO2 flux
    346              tr_seri_sum(:) = tr_seri(:,1,id_fco2_fos_fuel) + tr_seri(:,1,id_fco2_land_use) + &
    347                   tr_seri(:,1,id_fco2_land) + tr_seri(:,1,id_fco2_ocn)
    348              
    349              ! Calculate mean value of delta CO2 flux
    350              fco2_tmp(:) = tr_seri_sum(:) * airephy(:)
    351              CALL reduce_sum(SUM(fco2_tmp),sumtmp)
    352              delta_co2_ppm = sumtmp/airetot
    353              
    354              ! Add initial value for co2_ppm to delta value
    355              co2_ppm = delta_co2_ppm + co2_ppm0
    356           END IF
    357 
    358        ELSE IF (carbon_cycle_cpl) THEN ! no carbon_cycle_tr
    359 
    360           ! Calculate the total CO2 flux and integrate to get scalare value for the radiation scheme
    361           fco2_tmp(:) = (fos_fuel(:) + (fco2_lu_day(:) + fco2_land_day(:))*pctsrf(:,is_ter) &
    362                + fco2_ocn_day(:)*pctsrf(:,is_oce)) * fact
    363          
    364           ! Calculate mean value
    365           fco2_tmp(:) = fco2_tmp(:) * airephy(:)
    366           CALL reduce_sum(SUM(fco2_tmp),sumtmp)
    367           delta_co2_ppm = sumtmp/airetot
    368 
    369           ! Update current value of the atmospheric co2_ppm
    370           co2_ppm = co2_ppm + delta_co2_ppm
    371          
    372        END IF ! carbon_cycle_tr
    373 
    374        ! transformation of the atmospheric CO2 concentration for the radiation code
    375        RCO2 = co2_ppm * 1.0e-06  * 44.011/28.97
    376 
    377     END IF
    378 
    379     ! Calculate CO2 flux to send to ocean and land models : PISCES and ORCHIDEE         
    380     IF (endday .AND. carbon_cycle_cpl) THEN
    381        
    382        IF (carbon_cycle_tr) THEN
    383           IF (ntr_co2==1) THEN
    384 
    385              co2_send(:) = tr_seri(:,1,id_fco2_tot) + co2_ppm0
    386 
    387           ELSE ! ntr_co2 == 4
    388 
    389              co2_send(:) = tr_seri_sum(:) + co2_ppm0
    390 
    391           END IF
     452          ! Sum all co2 tracers to get the total delta CO2 flux at first model layer
     453          fco2_tmp(:) = 0.
     454          DO it = 1, ntr_co2
     455             fco2_tmp(1:klon) = fco2_tmp(1:klon) + tr_seri(1:klon,1,co2trac(it)%id)
     456          END DO
     457          co2_send(1:klon) = fco2_tmp(1:klon) + co2_ppm0
    392458       ELSE
    393459          ! Send a scalare value in 2D variable to ocean and land model (PISCES and ORCHIDEE)
    394           co2_send(:) = co2_ppm
     460          co2_send(1:klon) = co2_ppm
    395461       END IF
    396462

  • LMDZ5/trunk/libf/phylmd/change_srf_frac_mod.F90

    r996 r1454  
    99!
    1010! Change Surface Fractions
    11 !
     11! Author J Ghattas 2008
     12
    1213  SUBROUTINE change_srf_frac(itime, dtime, jour, &
    1314       pctsrf, alb1, alb2, tsurf, u10m, v10m, pbl_tke)
    … …  
    7677    END SELECT
    7778
    78     IF (is_modified) THEN
     79
    7980!****************************************************************************************
    8081! 2)
    … …  
    8485!
    8586!****************************************************************************************
     87    IF (is_modified) THEN
    8688 
    8789! Test and exit if a fraction is negative
    … …  
    150152       CALL pbl_surface_newfrac(itime, pctsrf, pctsrf_old, tsurf, alb1, alb2, u10m, v10m, pbl_tke)
    151153
     154    ELSE
     155       ! No modifcation should be done
     156       pctsrf(:,:) = pctsrf_old(:,:)
     157
    152158    END IF ! is_modified
    153159

  • LMDZ5/trunk/libf/phylmd/cpl_mod.F90

    r1403 r1454  
    100100  SUBROUTINE cpl_init(dtime, rlon, rlat)
    101101    USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, fco2_ocn_day
     102    USE surface_data
    102103
    103104    INCLUDE "dimensions.h"
    … …  
    270271    ENDIF    ! is_sequential
    271272   
     273
     274!*************************************************************************************
     275! compatibility test
     276!
     277!*************************************************************************************
     278    IF (carbon_cycle_cpl .AND. version_ocean=='opa8') THEN
     279       abort_message='carbon_cycle_cpl does not work with opa8'
     280       CALL abort_gcm(modname,abort_message,1)
     281    END IF
     282
    272283  END SUBROUTINE cpl_init
    273284 

  • LMDZ5/trunk/libf/phylmd/phyetat0.F

    r1403 r1454  
    2121      USE infotrac
    2222      USE traclmdz_mod,    ONLY : traclmdz_from_restart
    23       USE carbon_cycle_mod,ONLY : carbon_cycle_tr, carbon_cycle_cpl
     23      USE carbon_cycle_mod,ONLY :
     24     &     carbon_cycle_tr, carbon_cycle_cpl, co2_send
    2425
    2526      IMPLICIT none
    … …  
    133134
    134135       
    135 
    136          IF( clesphy0(1).NE.tab_cntrl( 5 ) )  THEN
    137              clesphy0(1)=tab_cntrl( 5 )
    138          ENDIF
    139 
    140          IF( clesphy0(2).NE.tab_cntrl( 6 ) )  THEN
    141              clesphy0(2)=tab_cntrl( 6 )
    142          ENDIF
    143 
    144          IF( clesphy0(3).NE.tab_cntrl( 7 ) )  THEN
    145              clesphy0(3)=tab_cntrl( 7 )
    146          ENDIF
    147 
    148          IF( clesphy0(4).NE.tab_cntrl( 8 ) )  THEN
    149              clesphy0(4)=tab_cntrl( 8 )
    150          ENDIF
    151 
    152          IF( clesphy0(5).NE.tab_cntrl( 9 ) )  THEN
    153              clesphy0(5)=tab_cntrl( 9 )
    154          ENDIF
    155 
    156          IF( clesphy0(6).NE.tab_cntrl( 10 ) )  THEN
    157              clesphy0(6)=tab_cntrl( 10 )
    158          ENDIF
    159 
    160          IF( clesphy0(7).NE.tab_cntrl( 11 ) )  THEN
    161              clesphy0(7)=tab_cntrl( 11 )
    162          ENDIF
    163 
    164          IF( clesphy0(8).NE.tab_cntrl( 12 ) )  THEN
    165              clesphy0(8)=tab_cntrl( 12 )
    166          ENDIF
    167 
     136      clesphy0(1)=tab_cntrl( 5 )
     137      clesphy0(2)=tab_cntrl( 6 )
     138      clesphy0(3)=tab_cntrl( 7 )
     139      clesphy0(4)=tab_cntrl( 8 )
     140      clesphy0(5)=tab_cntrl( 9 )
     141      clesphy0(6)=tab_cntrl( 10 )
     142      clesphy0(7)=tab_cntrl( 11 )
     143      clesphy0(8)=tab_cntrl( 12 )
    168144
    169145c
    … …  
    10781054
    10791055         END DO
    1080          
    10811056         CALL traclmdz_from_restart(trs)
     1057
     1058         IF (carbon_cycle_cpl) THEN
     1059            ALLOCATE(co2_send(klon), stat=ierr)
     1060            IF (ierr /= 0) CALL abort_gcm
     1061     &           ('phyetat0','pb allocation co2_send',1)
     1062            CALL get_field("co2_send",co2_send,found)
     1063            IF (.NOT. found) THEN
     1064               PRINT*,"phyetat0: Le champ <co2_send> est absent"
     1065               PRINT*,"Initialisation uniforme a co2_ppm=",co2_ppm
     1066               co2_send(:) = co2_ppm
     1067            END IF
     1068         END IF
    10821069      END IF
    10831070

  • LMDZ5/trunk/libf/phylmd/phyredem.F

    r1403 r1454  
    1515      USE infotrac
    1616      USE control_mod
    17 
     17      USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, co2_send
    1818
    1919      IMPLICIT none
    … …  
    337337            CALL put_field("trs_"//tname(iiq),"",trs(:,it))
    338338         END DO
     339         IF (carbon_cycle_cpl) THEN
     340            IF (.NOT. ALLOCATED(co2_send)) THEN
     341               ! This is the case of create_etat0_limit, ce0l
     342               ALLOCATE(co2_send(klon))
     343               co2_send(:) = co2_ppm0
     344            END IF
     345            CALL put_field("co2_send","co2_ppm for coupling",co2_send)
     346         END IF
    339347      END IF
    340348

  • LMDZ5/trunk/libf/phylmd/physiq.F

    r1428 r1454  
    12501250cym Attention pbase pas initialise dans concvl !!!!
    12511251          pbase=0
    1252           paire_ter(:)=0.   
    12531252cIM 180608
    12541253c         pmflxr=0.
    … …  
    33763375     I                   cdragh,coefh,u1,v1,ftsol,pctsrf,
    33773376     I                   frac_impa, frac_nucl,
    3378      I                   pphis,airephy,dtime,itap)
     3377     I                   pphis,airephy,dtime,itap,
     3378     I                   rlon,rlat,qx(:,:,ivap),da,phi,mp,upwd,dnwd)
    33793379
    33803380

  • LMDZ5/trunk/libf/phylmd/phytrac.F90

    r1421 r1454  
    6666!--------
    6767  REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
    68   REAL,DIMENSION(klon,klev),INTENT(IN)   :: u       !
    69   REAL,DIMENSION(klon,klev),INTENT(IN)   :: v       !
     68  REAL,DIMENSION(klon,klev),INTENT(IN)   :: u       ! variable not used
     69  REAL,DIMENSION(klon,klev),INTENT(IN)   :: v       ! variable not used
    7070  REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
    7171  REAL,DIMENSION(klon,klev),INTENT(IN)   :: rh      ! humidite relative
    … …  
    118118!--------------
    119119!
    120   REAL,DIMENSION(klon,klev),INTENT(IN) :: cdragh ! coeff drag pour T et Q
     120  REAL,DIMENSION(klon),INTENT(IN)      :: cdragh ! coeff drag pour T et Q
    121121  REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh  ! coeff melange CL (m**2/s)
    122122  REAL,DIMENSION(klon),INTENT(IN)      :: yu1    ! vents au premier niveau
    … …  
    213213     SELECT CASE(type_trac)
    214214     CASE('lmdz')
    215 !IM ajout t_seri, pplay, sh    CALL traclmdz_init(pctsrf, ftsol, tr_seri, aerosol, lessivage)
    216         CALL traclmdz_init(pctsrf, ftsol, tr_seri, t_seri, pplay, sh, aerosol, lessivage)
     215        CALL traclmdz_init(pctsrf, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
    217216     CASE('inca')
    218217        source(:,:)=0.

  • LMDZ5/trunk/libf/phylmd/read_map2D.F90

    r1279 r1454  
    1111
    1212! Input arguments
    13   CHARACTER(len=20), INTENT(IN) :: filename     ! name of file to read
    14   CHARACTER(len=20), INTENT(IN) :: varname      ! name of variable in file
     13  CHARACTER(len=*), INTENT(IN) :: filename     ! name of file to read
     14  CHARACTER(len=*), INTENT(IN) :: varname      ! name of variable in file
    1515  INTEGER, INTENT(IN)           :: timestep     ! actual timestep
    1616  LOGICAL, INTENT(IN)           :: inverse      ! TRUE if latitude needs to be inversed
    … …  
    2727  REAL, DIMENSION(nbp_lon,nbp_lat) :: var_glo2D_tmp ! 2D global
    2828  REAL, DIMENSION(klon_glo)        :: var_glo1D     ! 1D global
    29 
     29  INCLUDE "iniprint.h"
    3030
    3131! Read variable from file. Done by master process MPI and master thread OpenMP
    3232  IF (is_mpi_root .AND. is_omp_root) THEN
    33      ierr = NF90_OPEN (filename, NF90_NOWRITE, nid)
    34      IF (ierr /= NF90_NOERR) CALL abort_gcm(modname,'Problem in opening file '//filename,1)
     33     ierr = NF90_OPEN(trim(filename), NF90_NOWRITE, nid)
     34     IF (ierr /= NF90_NOERR) CALL write_err_mess('Problem in opening file')
    3535
    36      ierr = NF90_INQ_VARID(nid, varname, nvarid)
    37      IF (ierr /= NF90_NOERR) CALL abort_gcm(modname, 'The variable '//varname//' is absent in file',1)
     36     ierr = NF90_INQ_VARID(nid, trim(varname), nvarid)
     37     IF (ierr /= NF90_NOERR) CALL write_err_mess('The variable is absent in file')
    3838     
    3939     start=(/1,1,timestep/)
    4040     count=(/nbp_lon,nbp_lat,1/)
    4141     ierr = NF90_GET_VAR(nid, nvarid, var_glo2D,start,count)
    42      IF (ierr /= NF90_NOERR) CALL abort_gcm(modname, 'Problem in reading varaiable '//varname,1)
     42     IF (ierr /= NF90_NOERR) CALL write_err_mess('Problem in reading varaiable')
     43
     44     ierr = NF90_CLOSE(nid)
     45     IF (ierr /= NF90_NOERR) CALL write_err_mess('Problem in closing file')
    4346
    4447     ! Inverse latitude order
    … …  
    5356     CALL grid2Dto1D_glo(var_glo2D,var_glo1D)
    5457
     58     WRITE(lunout,*) 'in read_map2D, filename = ', trim(filename)
     59     WRITE(lunout,*) 'in read_map2D, varname  = ', trim(varname)
     60     WRITE(lunout,*) 'in read_map2D, timestep = ', timestep
    5561  ENDIF
    5662
    … …  
    5864  CALL scatter(var_glo1D, varout)
    5965
     66  CONTAINS
     67    SUBROUTINE write_err_mess(err_mess)
     68
     69      CHARACTER(len=*), INTENT(IN) :: err_mess
     70      INCLUDE "iniprint.h"
     71     
     72      WRITE(lunout,*) 'Error in read_map2D, filename = ', trim(filename)
     73      WRITE(lunout,*) 'Error in read_map2D, varname  = ', trim(varname)
     74      WRITE(lunout,*) 'Error in read_map2D, timestep = ', timestep
     75
     76      CALL abort_gcm(modname, err_mess, 1)
     77
     78    END SUBROUTINE write_err_mess
     79
    6080END SUBROUTINE read_map2D

  • LMDZ5/trunk/libf/phylmd/surf_land_orchidee_mod.F90

    r1279 r1454  
    4343    USE mod_surf_para
    4444    USE mod_synchro_omp
    45    
    46 USE carbon_cycle_mod, ONLY : carbon_cycle_cpl, fco2_land_inst, fco2_lu_inst
     45    USE carbon_cycle_mod, ONLY : carbon_cycle_cpl
    4746
    4847!   
    … …  
    138137    INTEGER                                   :: error
    139138    REAL, DIMENSION(klon)                     :: swdown_vrai
    140     REAL, DIMENSION(klon)                     :: fco2_land_comp  ! sur grille compresse
    141     REAL, DIMENSION(klon)                     :: fco2_lu_comp    ! sur grille compresse
    142139    CHARACTER (len = 20)                      :: modname = 'surf_land_orchidee'
    143140    CHARACTER (len = 80)                      :: abort_message
    … …  
    348345       ENDIF
    349346!
    350 ! Allocate variables needed for carbon_cycle_mod
     347! carbon_cycle_cpl not possible with this interface and version of ORHCHIDEE
    351348!
    352349       IF (carbon_cycle_cpl) THEN
    353           IF (.NOT. ALLOCATED(fco2_land_inst)) THEN
    354              ALLOCATE(fco2_land_inst(klon),stat=error)
    355              IF (error /= 0)  CALL abort_gcm(modname,'Pb in allocation fco2_land_inst',1)
    356              
    357              ALLOCATE(fco2_lu_inst(klon),stat=error)
    358              IF(error /=0) CALL abort_gcm(modname,'Pb in allocation fco2_lu_inst',1)
    359           END IF
     350          abort_message='carbon_cycle_cpl not yet possible with this interface of ORCHIDEE'
     351          CALL abort_gcm(modname,abort_message,1)
    360352       END IF
    361353       
    … …  
    464456
    465457   
    466 ! JG : TEMPORAIRE!!!! Les variables fco2_land_comp et fco2_lu_comp seront plus tard en sortie d'ORCHIDEE
    467 !      ici mis a valeur quelquonque pour test. Ces variables sont sur la grille compresse avec uniquement des points de terres
    468 
    469     fco2_land_comp(:) = 1.
    470     fco2_lu_comp(:)   = 10.
    471 
    472 ! Decompress variables for the module carbon_cycle_mod
    473     IF (carbon_cycle_cpl) THEN
    474        fco2_land_inst(:)=0.
    475        fco2_lu_inst(:)=0.
    476        
    477        DO igrid = 1, knon
    478           ireal = knindex(igrid)
    479           fco2_land_inst(ireal) = fco2_land_comp(igrid)
    480           fco2_lu_inst(ireal)   = fco2_lu_comp(igrid)
    481        END DO
    482     END IF
    483 
    484458  END SUBROUTINE surf_land_orchidee
    485459!

  • LMDZ5/trunk/libf/phylmd/surf_land_orchidee_noopenmp_mod.F90

    r1279 r1454  
    138138    INTEGER                                   :: ij, jj, igrid, ireal, index
    139139    INTEGER                                   :: error
     140    INTEGER, SAVE                             :: nb_fields_cpl ! number of fields for the climate-carbon coupling (between ATM and ORCHIDEE).
     141    REAL, SAVE, ALLOCATABLE, DIMENSION(:,:)   :: fields_cpl    ! Fluxes for the climate-carbon coupling
    140142    REAL, DIMENSION(klon)                     :: swdown_vrai
    141     REAL, DIMENSION(klon)                     :: fco2_land_comp  ! sur grille compresse
    142     REAL, DIMENSION(klon)                     :: fco2_lu_comp    ! sur grille compresse
    143143    CHARACTER (len = 20)                      :: modname = 'surf_land_orchidee'
    144144    CHARACTER (len = 80)                      :: abort_message
    … …  
    210210 
    211211    IF (debut) THEN
     212! Test de coherence
     213#ifndef ORCH_NEW
     214       ! Compilation avec orchidee nouvelle version necessaire avec carbon_cycle_cpl=y
     215       IF (carbon_cycle_cpl) THEN
     216          abort_message='You must define preprossing key ORCH_NEW when running carbon_cycle_cpl=y'
     217          CALL abort_gcm(modname,abort_message,1)
     218       END IF
     219#endif
    212220       ALLOCATE(ktindex(knon))
    213221       IF ( .NOT. ALLOCATED(albedo_keep)) THEN
    … …  
    342350!
    343351! Allocate variables needed for carbon_cycle_mod
    344 !
     352       IF ( carbon_cycle_cpl ) THEN
     353          nb_fields_cpl=2
     354       ELSE
     355          nb_fields_cpl=1
     356       END IF
     357
     358
    345359       IF (carbon_cycle_cpl) THEN
    346           IF (.NOT. ALLOCATED(fco2_land_inst)) THEN
    347              ALLOCATE(fco2_land_inst(klon),stat=error)
    348              IF (error /= 0)  CALL abort_gcm(modname,'Pb in allocation fco2_land_inst',1)
    349              
    350              ALLOCATE(fco2_lu_inst(klon),stat=error)
    351              IF(error /=0) CALL abort_gcm(modname,'Pb in allocation fco2_lu_inst',1)
    352           END IF
     360          ALLOCATE(fco2_land_inst(klon),stat=error)
     361          IF (error /= 0)  CALL abort_gcm(modname,'Pb in allocation fco2_land_inst',1)
     362         
     363          ALLOCATE(fco2_lu_inst(klon),stat=error)
     364          IF(error /=0) CALL abort_gcm(modname,'Pb in allocation fco2_lu_inst',1)
    353365       END IF
     366
     367       ALLOCATE(fields_cpl(klon,nb_fields_cpl), stat = error)
     368       IF (error /= 0) CALL abort_gcm(modname,'Pb in allocation fields_cpl',1)
    354369
    355370    ENDIF                          ! (fin debut)
    … …  
    406421               evap, fluxsens, fluxlat, coastalflow, riverflow, &
    407422               tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
    408                lon_scat, lat_scat, q2m, t2m)
     423               lon_scat, lat_scat, q2m, t2m &
     424#ifdef ORCH_NEW
     425               , nb_fields_cpl, fields_cpl)
     426#else
     427               )
     428#endif
    409429
    410430#else         
    411           ! Interface for ORCHIDEE version 1.9 or later(1.9.2, 1.9.3, 1.9.4) compiled in parallel mode(with preprocessing flag CPP_MPI)
     431          ! Interface for ORCHIDEE version 1.9 or later(1.9.2, 1.9.3, 1.9.4, 1.9.5) compiled in parallel mode(with preprocessing flag CPP_MPI)
    412432          CALL intersurf_main (itime+itau_phy-1, iim, jjm+1, offset, knon, ktindex, &
    413433               orch_comm, dtime, lrestart_read, lrestart_write, lalo, &
    … …  
    418438               evap(1:knon), fluxsens(1:knon), fluxlat(1:knon), coastalflow(1:knon), riverflow(1:knon), &
    419439               tsol_rad(1:knon), tsurf_new(1:knon), qsurf(1:knon), albedo_out(1:knon,:), emis_new(1:knon), z0_new(1:knon), &
    420                lon_scat, lat_scat, q2m, t2m)
     440               lon_scat, lat_scat, q2m, t2m &
     441#ifdef ORCH_NEW
     442               , nb_fields_cpl, fields_cpl(1:knon,:))
     443#else
     444               )
     445#endif
    421446#endif
    422447         
    … …  
    431456
    432457    IF (knon /=0) THEN
    433    
    434458#ifndef CPP_MPI
    435459       ! Interface for ORCHIDEE compiled in sequential mode(without preprocessing flag CPP_MPI)
    … …  
    442466            evap, fluxsens, fluxlat, coastalflow, riverflow, &
    443467            tsol_rad, tsurf_new, qsurf, albedo_out, emis_new, z0_new, &
    444             lon_scat, lat_scat, q2m, t2m)
    445        
     468            lon_scat, lat_scat, q2m, t2m &
     469#ifdef ORCH_NEW
     470            , nb_fields_cpl, fields_cpl)
     471#else
     472            )
     473#endif
    446474#else
    447475       ! Interface for ORCHIDEE version 1.9 or later compiled in parallel mode(with preprocessing flag CPP_MPI)
    … …  
    454482            evap(1:knon), fluxsens(1:knon), fluxlat(1:knon), coastalflow(1:knon), riverflow(1:knon), &
    455483            tsol_rad(1:knon), tsurf_new(1:knon), qsurf(1:knon), albedo_out(1:knon,:), emis_new(1:knon), z0_new(1:knon), &
    456             lon_scat, lat_scat, q2m, t2m)
    457 #endif
    458        
     484            lon_scat, lat_scat, q2m, t2m &
     485#ifdef ORCH_NEW
     486            , nb_fields_cpl, fields_cpl(1:knon,:))
     487#else
     488            )
     489#endif
     490#endif
    459491    ENDIF
    460492
    … …  
    478510
    479511    IF (debut) lrestart_read = .FALSE.
    480 
    481 
    482 ! JG : TEMPORAIRE!!!! Les variables fco2_land_comp et fco2_lu_comp seront plus tard en sortie d'ORCHIDEE
    483 !      ici mis a valeur quelquonque pour test. Ces variables sont sur la grille compresse avec uniquement des points de terres
    484 
    485     fco2_land_comp(:) = 1.
    486     fco2_lu_comp(:)   = 10.
    487512
    488513! Decompress variables for the module carbon_cycle_mod
    … …  
    493518       DO igrid = 1, knon
    494519          ireal = knindex(igrid)
    495           fco2_land_inst(ireal) = fco2_land_comp(igrid)
    496           fco2_lu_inst(ireal)   = fco2_lu_comp(igrid)
     520          fco2_land_inst(ireal) = fields_cpl(igrid,1)
     521          fco2_lu_inst(ireal)   = fields_cpl(igrid,2)
    497522       END DO
    498523    END IF

  • LMDZ5/trunk/libf/phylmd/traclmdz_mod.F90

    r1421 r1454  
    8484
    8585
    86   SUBROUTINE traclmdz_init(pctsrf, ftsol, tr_seri, t_seri, pplay, sh, aerosol, lessivage)
     86  SUBROUTINE traclmdz_init(pctsrf, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
    8787    ! This subroutine allocates and initialize module variables and control variables.
    8888    ! Initialization of the tracers should be done here only for those not found in the restart file.
    … …  
    104104    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay   ! pression pour le mileu de chaque couche (en Pa)
    105105    REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
     106    REAL,INTENT(IN)                        :: pdtphys ! Pas d'integration pour la physique (seconde) 
    106107
    107108! Output variables
    … …  
    226227! ----------------------------------------------
    227228    IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
    228        CALL carbon_cycle_init(tr_seri, aerosol, radio)
     229       CALL carbon_cycle_init(tr_seri, pdtphys, aerosol, radio)
    229230    END IF
    230231
    … …  
    312313!--------------
    313314!
    314     REAL,DIMENSION(klon,klev),INTENT(IN) :: cdragh     ! coeff drag pour T et Q
     315    REAL,DIMENSION(klon),INTENT(IN)      :: cdragh     ! coeff drag pour T et Q
    315316    REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh      ! coeff melange CL (m**2/s)
    316317    REAL,DIMENSION(klon),INTENT(IN)      :: yu1        ! vents au premier niveau
    … …  
    546547!======================================================================
    547548    IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
    548        CALL carbon_cycle(nstep, pdtphys, pctsrf, tr_seri)
     549       CALL carbon_cycle(nstep, pdtphys, pctsrf, tr_seri, source)
    549550    END IF
    550551
Note: See TracChangeset for help on using the changeset viewer.