Ignore:
Timestamp:
Mar 3, 2015, 2:41:13 PM (10 years ago)
Author:
Laurent Fairhead
Message:

Merged trunk changes -r2186:2216 into testing branch

Location:
LMDZ5/branches/testing
Files:
69 edited
3 copied

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/testing

  • LMDZ5/branches/testing/libf/bibio/i1mach.F

    r1910 r2220  
    11*DECK I1MACH
    22      INTEGER FUNCTION I1MACH (I)
     3      IMPLICIT NONE
    34C***BEGIN PROLOGUE  I1MACH
    45C***PURPOSE  Return integer machine dependent constants.
    … …  
    9596      SAVE IMACH
    9697      EQUIVALENCE (IMACH(4),OUTPUT)
     98      INTEGER I
    9799C***FIRST EXECUTABLE STATEMENT  I1MACH
    98100      IMACH( 1) =         5

  • LMDZ5/branches/testing/libf/bibio/j4save.F

    r1910 r2220  
    11*DECK J4SAVE
    22      FUNCTION J4SAVE (IWHICH, IVALUE, ISET)
     3      IMPLICIT NONE
    34C***BEGIN PROLOGUE  J4SAVE
    45C***SUBSIDIARY
    … …  
    5960      DATA IPARAM(5)/1/
    6061      DATA IPARAM(6),IPARAM(7),IPARAM(8),IPARAM(9)/0,0,0,0/
     62      INTEGER J4SAVE,IWHICH,IVALUE
    6163C***FIRST EXECUTABLE STATEMENT  J4SAVE
    6264      J4SAVE = IPARAM(IWHICH)

  • LMDZ5/branches/testing/libf/bibio/xercnt.F

    r1910 r2220  
    11*DECK XERCNT
    22      SUBROUTINE XERCNT (LIBRAR, SUBROU, MESSG, NERR, LEVEL, KONTRL)
     3      IMPLICIT NONE
    34C***BEGIN PROLOGUE  XERCNT
    45C***SUBSIDIARY
    … …  
    5657C***END PROLOGUE  XERCNT
    5758      CHARACTER*(*) LIBRAR, SUBROU, MESSG
     59      INTEGER NERR, LEVEL, KONTRL
    5860C***FIRST EXECUTABLE STATEMENT  XERCNT
    5961      RETURN

  • LMDZ5/branches/testing/libf/bibio/xermsg.F

    r1910 r2220  
    11*DECK XERMSG
    22      SUBROUTINE XERMSG (LIBRAR, SUBROU, MESSG, NERR, LEVEL)
     3      IMPLICIT NONE
    34C***BEGIN PROLOGUE  XERMSG
    45C***PURPOSE  Process error messages for SLATEC and other libraries.
    … …  
    189190      CHARACTER*72  TEMP
    190191      CHARACTER*20  LFIRST
     192      INTEGER NERR, LEVEL, LKNTRL
     193      INTEGER J4SAVE, MAXMES, KDUMMY, I, KOUNT, LERR, LLEVEL
     194      INTEGER MKNTRL, LTEMP
    191195C***FIRST EXECUTABLE STATEMENT  XERMSG
    192196      LKNTRL = J4SAVE (2, 0, .FALSE.)

  • LMDZ5/branches/testing/libf/bibio/xerprn.F

    r1910 r2220  
    11*DECK XERPRN
    22      SUBROUTINE XERPRN (PREFIX, NPREF, MESSG, NWRAP)
     3      IMPLICIT NONE
    34C***BEGIN PROLOGUE  XERPRN
    45C***SUBSIDIARY
    … …  
    8182      CHARACTER*2 NEWLIN
    8283      PARAMETER (NEWLIN = '$$')
     84      INTEGER N, I1MACH, I, LPREF, LWRAP, LENMSG, NEXTC
     85      INTEGER LPIECE, IDELTA
    8386C***FIRST EXECUTABLE STATEMENT  XERPRN
    8487      CALL XGETUA(IU,NUNIT)

  • LMDZ5/branches/testing/libf/bibio/xersve.F

    r1910 r2220  
    22      SUBROUTINE XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL,
    33     +   ICOUNT)
     4      IMPLICIT NONE
    45C***BEGIN PROLOGUE  XERSVE
    56C***SUBSIDIARY
    … …  
    5859C   920501  Reformatted the REFERENCES section.  (WRB)
    5960C***END PROLOGUE  XERSVE
    60       PARAMETER (LENTAB=10)
     61      INTEGER,PARAMETER :: LENTAB=10
    6162      INTEGER LUN(5)
    6263      CHARACTER*(*) LIBRAR, SUBROU, MESSG
    … …  
    6667      SAVE LIBTAB, SUBTAB, MESTAB, NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG
    6768      DATA KOUNTX/0/, NMSG/0/
     69      INTEGER NERR,LEVEL,KONTRL
     70      INTEGER NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG
     71      INTEGER KFLAG, ICOUNT, NUNIT, KUNIT, IUNIT, I1MACH, I
    6872C***FIRST EXECUTABLE STATEMENT  XERSVE
    6973C

  • LMDZ5/branches/testing/libf/bibio/xgetua.F

    r1910 r2220  
    11*DECK XGETUA
    22      SUBROUTINE XGETUA (IUNITA, N)
     3      IMPLICIT NONE
    34C***BEGIN PROLOGUE  XGETUA
    45C***PURPOSE  Return unit number(s) to which error messages are being
    … …  
    4142C***END PROLOGUE  XGETUA
    4243      DIMENSION IUNITA(5)
     44      INTEGER IUNITA, N, J4SAVE, INDEX, I
    4345C***FIRST EXECUTABLE STATEMENT  XGETUA
    4446      N = J4SAVE(5,0,.FALSE.)

  • LMDZ5/branches/testing/libf/dyn3d_common/grid_atob.F

    r2160 r2220  
    936936c
    937937      SUBROUTINE dist_sphe(rf_lon,rf_lat,rlon,rlat,im,jm,distance)
     938      IMPLICIT NONE
    938939c
    939940c Auteur: Laurent Li (le 30 decembre 1996)
    … …  
    958959      REAL radius
    959960      PARAMETER (radius=6371229.)
     961      INTEGER i,j
    960962c
    961963      pi = 4.0 * ATAN(1.0)

  • LMDZ5/branches/testing/libf/dyn3d_common/grid_noro.F

    r1999 r2220  
    5454C=======================================================================
    5555
    56       IMPLICIT INTEGER (I,J)
    57       IMPLICIT REAL(X,Z)
     56      IMPLICIT NONE
    5857     
    59           parameter(iusn=2160,jusn=1080,iext=216, epsfra = 1.e-5)
     58      INTEGER,PARAMETER :: iusn=2160,jusn=1080,iext=216
     59      REAL,PARAMETER :: epsfra = 1.e-5
    6060#include "dimensions.h"
    6161          REAL xusn(iusn+2*iext),yusn(jusn+2)   
    … …  
    8989      REAL a(2200),b(2200),c(1100),d(1100)
    9090      logical masque_lu
     91      INTEGER i,j,ii,jj
     92      REAL xpi, rad, zdeltay, masque
     93      REAL zdeltax, zlenx, zleny, xincr
     94      REAL zbordnor, zbordsud, weighy, zbordest, zbordoue, weighx
     95      REAL zllmmea, zllmstd, zllmsig, zllmgam, zllmpic, zllmval
     96      REAL zllmthe, zminthe
     97      REAL xk, xl, xm, xp, xq, xw
     98      REAL zmeanor, zmeasud, zstdnor, zstdsud, zsignor, zsigsud
     99      REAL zweinor, zweisud, zpicnor, zpicsud, zvalnor, zvalsud
    91100c
    92101      print *,' parametres de l orographie a l echelle sous maille'
    … …  
    455464
    456465      SUBROUTINE MVA9(X,IMAR,JMAR)
    457 
     466      IMPLICIT NONE
    458467C MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
    459 
     468      INTEGER IMAR,JMAR
    460469      REAL X(IMAR,JMAR),XF(IMAR,JMAR)
    461470      real WEIGHTpb(-1:1,-1:1)
     471      INTEGER I,J,IS,JS
     472      REAL SUM
    462473
    463474

  • LMDZ5/branches/testing/libf/dyn3d_common/juldate.F

    r1999 r2220  
    77c       En entree:an,mois,jour,heure,min.,sec.
    88c       En sortie:tjd
    9         implicit real (a-h,o-z)
     9        IMPLICIT NONE
     10        INTEGER,INTENT(IN) :: ian,imoi,ijou,oh,om,os
     11        REAL,INTENT(OUT) :: tjd,tjdsec
     12       
     13        REAL frac,year,rmon,cf,a,b
     14        INTEGER ojou
     15       
    1016        frac=((os/60.+om)/60.+oh)/24.
    1117        ojou=dble(ijou)+frac

  • LMDZ5/branches/testing/libf/dyn3d_common/massbar.F

    r1999 r2220  
    33!
    44      SUBROUTINE massbar(  masse, massebx, masseby )
     5      IMPLICIT NONE
    56c
    67c **********************************************************************
    … …  
    2425      REAL    masse( ip1jmp1,llm ), massebx( ip1jmp1,llm )  ,
    2526     *      masseby(   ip1jm,llm )
     27      INTEGER ij,l
    2628c
    2729c

  • LMDZ5/branches/testing/libf/dyn3d_common/massbarxy.F

    r1999 r2220  
    33!
    44      SUBROUTINE massbarxy(  masse, massebxy )
     5      IMPLICIT NONE
    56c
    67c **********************************************************************
    … …  
    2324c
    2425       REAL  masse( ip1jmp1,llm ), massebxy( ip1jm,llm )
     26       INTEGER ij,l
    2527c
    2628

  • LMDZ5/branches/testing/libf/dyn3d_common/ppm3d.F

    r1999 r2220  
    6666     &                  JNP,j1,NLAY,AP,BP,PT,AE,fill,dum,Umax)
    6767
    68 c      implicit none
     68      implicit none
    6969
    7070c     rajout de déclarations
    … …  
    270270C     User modifiable parameters
    271271C
    272       parameter (Jmax = 361, kmax = 150)
     272      integer,parameter :: Jmax = 361, kmax = 150
    273273C
    274274C ****6***0*********0*********0*********0*********0*********0**********72
    … …  
    299299      data NDT0, NSTEP /0, 0/
    300300      data cross /.true./
     301      REAL DTDY, DTDY5, RCAP
     302      INTEGER JS0, JN0, IML, JMR, IMJM
    301303      SAVE DTDY, DTDY5, RCAP, JS0, JN0, IML,
    302304     &     DTDX, DTDX5, ACOSP, COSP, COSE, DAP,DBK
    303305C
     306      INTEGER NDT0, NSTEP, j2, k,j,i,ic,l,JS,JN,IMH
     307      INTEGER IU,IIU,JT,iad,jad,krd
     308      REAL r23,r3,PI,DL,DP,DT,CR1,MAXDT,ZTC,D5
     309      REAL sum1,sum2,ru
    304310           
    305311      JMR = JNP -1
    … …  
    756762      subroutine FZPPM(IMR,JNP,NLAY,j1,DQ,WZ,P,DC,DQDT,AR,AL,A6,
    757763     &                 flux,wk1,wk2,wz2,delp,KORD)
    758       parameter ( kmax = 150 )
    759       parameter ( R23 = 2./3., R3 = 1./3.)
     764      implicit none
     765      integer,parameter :: kmax = 150
     766      real,parameter :: R23 = 2./3., R3 = 1./3.
     767      integer IMR,JNP,NLAY,J1,KORD
    760768      real WZ(IMR,JNP,NLAY),P(IMR,JNP,NLAY),DC(IMR,JNP,NLAY),
    761769     &     wk1(IMR,*),delp(IMR,JNP,NLAY),DQ(IMR,JNP,NLAY),
    … …  
    764772      real AR(IMR,*),AL(IMR,*),A6(IMR,*),flux(IMR,*),wk2(IMR,*),
    765773     &     wz2(IMR,*)
     774      integer JMR,IMJM,NLAYM1,LMT,K,I,J
     775      real c0,c1,c2,tmp,qmax,qmin,a,b,fct,a1,a2,cm,cp
    766776C
    767777      JMR = JNP - 1
    … …  
    922932      subroutine xtp(IMR,JNP,IML,j1,j2,JN,JS,PU,DQ,Q,UC,
    923933     &               fx1,xmass,IORD)
     934      implicit none
     935      integer IMR,JNP,IML,j1,j2,JN,JS,IORD
     936      real PU,DQ,Q,UC,fx1,xmass
     937      real dc,qtmp
     938      integer ISAVE(IMR)
    924939      dimension UC(IMR,*),DC(-IML:IMR+IML+1),xmass(IMR,JNP)
    925940     &    ,fx1(IMR+1),DQ(IMR,JNP),qtmp(-IML:IMR+1+IML)
    926       dimension PU(IMR,JNP),Q(IMR,JNP),ISAVE(IMR)
     941      dimension PU(IMR,JNP),Q(IMR,JNP)
     942      integer jvan,j1vl,j2vl,j,i,iu,itmp,ist,imp
     943      real rut
    927944C
    928945      IMP = IMR + 1
    … …  
    10311048C
    10321049      subroutine fxppm(IMR,IML,UT,P,DC,flux,IORD)
    1033       parameter ( R3 = 1./3., R23 = 2./3. )
     1050      implicit none
     1051      integer IMR,IML,IORD
     1052      real UT,P,DC,flux
     1053      real,parameter ::  R3 = 1./3., R23 = 2./3.
    10341054      DIMENSION UT(*),flux(*),P(-IML:IMR+IML+1),DC(-IML:IMR+IML+1)
    1035       DIMENSION AR(0:IMR),AL(0:IMR),A6(0:IMR)
    1036       integer LMT 
     1055      REAL :: AR(0:IMR),AL(0:IMR),A6(0:IMR)
     1056      integer LMT,IMP,JLVL,i
    10371057c      logical first
    10381058c      data first /.true./
    … …  
    10881108C
    10891109      subroutine xmist(IMR,IML,P,DC)
    1090       parameter( R24 = 1./24.)
    1091       dimension P(-IML:IMR+1+IML),DC(-IML:IMR+1+IML)
     1110      implicit none
     1111      integer IMR,IML
     1112      real,parameter :: R24 = 1./24.
     1113      real :: P(-IML:IMR+1+IML),DC(-IML:IMR+1+IML)
     1114      integer :: i
     1115      real :: tmp,pmax,pmin
    10921116C
    10931117      do 10  i=1,IMR
    … …  
    11011125      subroutine ytp(IMR,JNP,j1,j2,acosp,RCAP,DQ,P,VC,DC2
    11021126     &              ,ymass,fx,A6,AR,AL,JORD)
     1127      implicit none
     1128      integer :: IMR,JNP,j1,j2,JORD
     1129      real :: acosp,RCAP,DQ,P,VC,DC2,ymass,fx,A6,AR,AL
    11031130      dimension P(IMR,JNP),VC(IMR,JNP),ymass(IMR,JNP)
    11041131     &       ,DC2(IMR,JNP),DQ(IMR,JNP),acosp(JNP)
    11051132C Work array
    11061133      DIMENSION fx(IMR,JNP),AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP)
     1134      integer :: JMR,len,i,jt,j
     1135      real :: sum1,sum2
    11071136C
    11081137      JMR = JNP - 1
    … …  
    11611190C
    11621191      subroutine  ymist(IMR,JNP,j1,P,DC,ID)
    1163       parameter ( R24 = 1./24. )
    1164       dimension P(IMR,JNP),DC(IMR,JNP)
     1192      implicit none
     1193      integer :: IMR,JNP,j1,ID
     1194      real,parameter :: R24 = 1./24.
     1195      real :: P(IMR,JNP),DC(IMR,JNP)
     1196      integer :: iimh,jmr,ijm3,imh,i
     1197      real :: pmax,pmin,tmp
    11651198C
    11661199      IMH = IMR / 2
    … …  
    12391272C
    12401273      subroutine fyppm(VC,P,DC,flux,IMR,JNP,j1,j2,A6,AR,AL,JORD)
    1241       parameter ( R3 = 1./3., R23 = 2./3. )
     1274      implicit none
     1275      integer IMR,JNP,j1,j2,JORD
     1276      real,parameter :: R3 = 1./3., R23 = 2./3.
    12421277      real VC(IMR,*),flux(IMR,*),P(IMR,*),DC(IMR,*)
    12431278C Local work arrays.
    12441279      real AR(IMR,JNP),AL(IMR,JNP),A6(IMR,JNP)
    1245       integer LMT
     1280      integer LMT,i
     1281      integer IMH,JMR,j11,IMJM1,len
    12461282c      logical first
    12471283C      data first /.true./
    … …  
    13151351C
    13161352        subroutine yadv(IMR,JNP,j1,j2,p,VA,ady,wk,IAD)
     1353        implicit none
     1354        integer IMR,JNP,j1,j2,IAD
    13171355        REAL p(IMR,JNP),ady(IMR,JNP),VA(IMR,JNP)
    13181356        REAL WK(IMR,-1:JNP+2)
     1357        INTEGER JMR,IMH,i,j,jp
     1358        REAL rv,a1,b1,sum1,sum2
    13191359C
    13201360        JMR = JNP-1
    … …  
    14011441C
    14021442        subroutine xadv(IMR,JNP,j1,j2,p,UA,JS,JN,IML,adx,IAD)
     1443        implicit none
     1444        INTEGER IMR,JNP,j1,j2,JS,JN,IML,IAD
    14031445        REAL p(IMR,JNP),adx(IMR,JNP),qtmp(-IMR:IMR+IMR),UA(IMR,JNP)
     1446        INTEGER JMR,j,i,ip,iu,iiu
     1447        REAL ru,a1,b1
    14041448C
    14051449        JMR = JNP-1
    … …  
    14891533C
    14901534      subroutine lmtppm(DC,A6,AR,AL,P,IM,LMT)
     1535      implicit none
    14911536C
    14921537C A6 =  CURVATURE OF THE TEST PARABOLA
    … …  
    15031548C LMT = 2: POSITIVE-DEFINITE CONSTRAINT
    15041549C
    1505       parameter ( R12 = 1./12. )
    1506       dimension A6(IM),AR(IM),AL(IM),P(IM),DC(IM)
     1550      real,parameter :: R12 = 1./12.
     1551      real :: A6(IM),AR(IM),AL(IM),P(IM),DC(IM)
     1552      integer :: IM,LMT
     1553      INTEGER i
     1554      REAL da1,da2,a6da,fmin
    15071555C
    15081556      if(LMT.eq.0) then
    … …  
    15641612C
    15651613      subroutine A2C(U,V,IMR,JMR,j1,j2,CRX,CRY,dtdx5,DTDY5)
    1566       dimension U(IMR,*),V(IMR,*),CRX(IMR,*),CRY(IMR,*),DTDX5(*)
     1614      implicit none
     1615      integer IMR,JMR,j1,j2
     1616      real :: U(IMR,*),V(IMR,*),CRX(IMR,*),CRY(IMR,*),DTDX5(*),DTDY5
     1617      integer i,j
    15671618C
    15681619      do 35 j=j1,j2
    … …  
    15791630C
    15801631      subroutine cosa(cosp,cose,JNP,PI,DP)
    1581       dimension cosp(*),cose(*)
     1632      implicit none
     1633      integer JNP
     1634      real :: cosp(*),cose(*),PI,DP
     1635      integer JMR,j,jeq
     1636      real ph5
    15821637      JMR = JNP-1
    15831638      do 55 j=2,JNP
    … …  
    16061661C
    16071662      subroutine cosc(cosp,cose,JNP,PI,DP)
    1608       dimension cosp(*),cose(*)
     1663      implicit none
     1664      integer JNP
     1665      real :: cosp(*),cose(*),PI,DP
     1666      real phi
     1667      integer j
    16091668C
    16101669      phi = -0.5*PI
    … …  
    16281687     &                   cross,IC,NSTEP)
    16291688C
    1630       parameter( tiny = 1.E-60 )
    1631       DIMENSION Q(IMR,JNP,NLAY),qtmp(IMR,JNP),cosp(*),acosp(*)
     1689      real,parameter :: tiny = 1.E-60
     1690      INTEGER :: IMR,JNP,NLAY,j1,j2,IC,NSTEP
     1691      REAL :: Q(IMR,JNP,NLAY),qtmp(IMR,JNP),cosp(*),acosp(*)
    16321692      logical cross
     1693      INTEGER :: NLAYM1,len,ip,L,icr,ipy,ipx,i
     1694      real :: qup,qly,dup,sum
    16331695C
    16341696      NLAYM1 = NLAY-1
    … …  
    17301792C
    17311793      subroutine filcr(q,IMR,JNP,j1,j2,cosp,acosp,icr,tiny)
    1732       dimension q(IMR,*),cosp(*),acosp(*)
     1794      implicit none
     1795      integer :: IMR,JNP,j1,j2,icr
     1796      real :: q(IMR,*),cosp(*),acosp(*),tiny
     1797      integer :: i,j
     1798      real :: dq,dn,d0,d1,ds,d2
    17331799      icr = 0
    17341800      do 65 j=j1+1,j2-1
    … …  
    18281894C
    18291895      subroutine filns(q,IMR,JNP,j1,j2,cosp,acosp,ipy,tiny)
    1830       dimension q(IMR,*),cosp(*),acosp(*)
     1896      implicit none
     1897      integer :: IMR,JNP,j1,j2,ipy
     1898      real :: q(IMR,*),cosp(*),acosp(*),tiny
     1899      real :: DP,CAP1,dq,dn,d0,d1,ds,d2
     1900      INTEGER :: i,j
    18311901c      logical first
    18321902c      data first /.true./
    … …  
    19101980C
    19111981      subroutine filew(q,qtmp,IMR,JNP,j1,j2,ipx,tiny)
    1912       dimension q(IMR,*),qtmp(JNP,IMR)
     1982      implicit none
     1983      integer :: IMR,JNP,j1,j2,ipx
     1984      real :: q(IMR,*),qtmp(JNP,IMR),tiny
     1985      integer :: i,j
     1986      real :: d0,d1,d2
    19131987C
    19141988      ipx = 0
    … …  
    19832057C
    19842058      subroutine zflip(q,im,km,nc)
     2059      implicit none
    19852060C This routine flip the array q (in the vertical).
     2061      integer :: im,km,nc
    19862062      real q(im,km,nc)
    19872063C local dynamic array
    19882064      real qtmp(im,km)
     2065      integer IC,k,i
    19892066C
    19902067      do 4000 IC = 1, nc

  • LMDZ5/branches/testing/libf/dyn3d_common/ran1.F

    r1999 r2220  
    33!
    44      FUNCTION RAN1(IDUM)
    5       DIMENSION R(97)
    6       save r
    7       save iff,ix1,ix2,ix3
    8       PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.8580247E-6)
    9       PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.4373773E-6)
    10       PARAMETER (M3=243000,IA3=4561,IC3=51349)
    11       DATA IFF /0/
     5      IMPLICIT NONE
     6      REAL RAN1
     7      REAL,SAVE :: R(97)
     8      REAL,PARAMETER :: RM1=3.8580247E-6,RM2=7.4373773E-6
     9      INTEGER,SAVE :: IFF=0
     10      integer,save :: ix1,ix2,ix3
     11      INTEGER,PARAMETER :: M1=259200,IA1=7141,IC1=54773
     12      INTEGER,PARAMETER :: M2=134456,IA2=8121,IC2=28411
     13      INTEGER,PARAMETER :: M3=243000,IA3=4561,IC3=51349
     14      INTEGER :: IDUM,J
     15
    1216      IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
    1317        IFF=1

  • LMDZ5/branches/testing/libf/filtrez/acc.F

    r1910 r2220  
    33!
    44        subroutine acc(vec,d,im)
    5         dimension vec(im,im),d(im)
     5        implicit none
     6        integer :: im
     7        real :: vec(im,im),d(im)
     8        integer :: i,j
     9        real ::sum
     10        real,external :: ssum
    611        do j=1,im
    712          do i=1,im

  • LMDZ5/branches/testing/libf/filtrez/eigen.F

    r1910 r2220  
    33!
    44      SUBROUTINE eigen( e,d)
     5      IMPLICIT NONE
    56#include "dimensions.h"
    6       dimension e( iim,iim ), d( iim )
    7       dimension asm( iim )
     7      real :: e( iim,iim ), d( iim )
     8      real :: asm( iim )
     9      integer :: im,i,j
    810      im=iim
    911c

  • LMDZ5/branches/testing/libf/grid/dimension/makdim

    r1910 r2220  
    55if (( $# < 1 )) || (( $# > 3 ))
    66then
    7  echo "Wrong number of parameters in $0 !!!"
    8  echo " Usage:"
    9  echo "  $0 [im] [jm] lm"
    10  echo " where im, jm and lm are the dimensions"
    11  exit
     7    echo "Wrong number of parameters in $0 !!!"
     8    echo " Usage:"
     9    echo "  $0 [im] [jm] lm"
     10    echo " where im, jm and lm are the dimensions"
     11    exit 1
     12fi
     13
     14if (($1 % 8 != 0)) && (( $# = 3 ))
     15then
     16    echo "The number of longitudes must be a multiple of 8."
     17    echo "See the files dyn3d/groupe.F and dyn3dmem/groupe_loc.F."
     18    exit 1
    1219fi
    1320
    1421# build "fichnom", the relevant 'dimensions.im.jm.lm' file name
    1522for i in $*
    16   do
    17   list=$list.$i
     23do
     24    list=$list.$i
    1825done
    1926fichdim=dimensions${list}
    2027
    2128if [ ! -f $fichdim ]
    22     then
     29then
    2330#    echo "$fichdim does not exist"
    2431
    2532    # assign values of im, jm and lm
    2633    if [ $# -ge 3 ]
    27         then
     34    then
    2835        im=$1
    2936        jm=$2
    … …  
    3138        ndm=1
    3239    elif [ $# -ge 2 ]
    33         then
     40    then
    3441        im=1
    3542        jm=$1
    … …  
    3744        ndm=1
    3845    elif [ $# -ge 1 ]
    39         then
     46    then
    4047        im=1
    4148        jm=1
    … …  
    4552
    4653# since the file doesn't exist, we create it
    47 cat << EOF > $fichdim
     54    cat << EOF > $fichdim
    4855!-----------------------------------------------------------------------
    4956!   INCLUDE 'dimensions.h'
    … …  
    6471# remove 'old' dimensions.h file (if any) and replace it with new one
    6572if [ -f ../dimensions.h ] ; then
    66 \rm ../dimensions.h
     73    \rm ../dimensions.h
    6774fi
    6875tar cf - $fichdim | ( cd .. ; tar xf - ; mv $fichdim dimensions.h )

  • LMDZ5/branches/testing/libf/phylmd/1DUTILS.h

    r2187 r2220  
    9999!             LS convergence imposed from  RICO (cst)
    100100!         = 6 ==> forcing_amma = .true.
     101!         = 10 ==> forcing_case = .true.
     102!             initial profiles from case.nc file
    101103!         = 40 ==> forcing_GCSSold = .true.
    102104!             initial profile from GCSS file
    … …  
    132134          CALL getin('turb_fcg',xTurb_fcg_gcssold)
    133135        ENDIF
     136
     137!Paramètres de forçage
     138!Config  Key  = tend_t
     139!Config  Desc = forcage ou non par advection de T
     140!Config  Def  = false
     141!Config  Help = forcage ou non par advection de T
     142       tend_t =0
     143       CALL getin('tend_t',tend_t)
     144
     145!Config  Key  = tend_q
     146!Config  Desc = forcage ou non par advection de q
     147!Config  Def  = false
     148!Config  Help = forcage ou non par advection de q
     149       tend_q =0
     150       CALL getin('tend_q',tend_q)
     151
     152!Config  Key  = tend_u
     153!Config  Desc = forcage ou non par advection de u
     154!Config  Def  = false
     155!Config  Help = forcage ou non par advection de u
     156       tend_u =0
     157       CALL getin('tend_u',tend_u)
     158
     159!Config  Key  = tend_v
     160!Config  Desc = forcage ou non par advection de v
     161!Config  Def  = false
     162!Config  Help = forcage ou non par advection de v
     163       tend_v =0
     164       CALL getin('tend_v',tend_v)
     165
     166!Config  Key  = tend_w
     167!Config  Desc = forcage ou non par vitesse verticale
     168!Config  Def  = false
     169!Config  Help = forcage ou non par vitesse verticale
     170       tend_w =0
     171       CALL getin('tend_w',tend_w)
     172
     173!Config  Key  = tend_rayo
     174!Config  Desc = forcage ou non par dtrad
     175!Config  Def  = false
     176!Config  Help = forcage ou non par dtrad
     177       tend_rayo =0
     178       CALL getin('tend_rayo',tend_rayo)
     179
     180
     181!Config  Key  = nudge_t
     182!Config  Desc = constante de nudging de T
     183!Config  Def  = false
     184!Config  Help = constante de nudging de T
     185       nudge_t =0.
     186       CALL getin('nudge_t',nudge_t)
     187
     188!Config  Key  = nudge_q
     189!Config  Desc = constante de nudging de q
     190!Config  Def  = false
     191!Config  Help = constante de nudging de q
     192       nudge_q =0.
     193       CALL getin('nudge_q',nudge_q)
     194
     195!Config  Key  = nudge_u
     196!Config  Desc = constante de nudging de u
     197!Config  Def  = false
     198!Config  Help = constante de nudging de u
     199       nudge_u =0.
     200       CALL getin('nudge_u',nudge_u)
     201
     202!Config  Key  = nudge_v
     203!Config  Desc = constante de nudging de v
     204!Config  Def  = false
     205!Config  Help = constante de nudging de v
     206       nudge_v =0.
     207       CALL getin('nudge_v',nudge_v)
     208
     209!Config  Key  = nudge_w
     210!Config  Desc = constante de nudging de w
     211!Config  Def  = false
     212!Config  Help = constante de nudging de w
     213       nudge_w =0.
     214       CALL getin('nudge_w',nudge_w)
     215
    134216
    135217!Config  Key  = iflag_nudge
    … …  
    24312513 
    24322514!=====================================================================
     2515       SUBROUTINE interp_case_vertical(play,nlev_cas,plev_prof_cas            &
     2516     &         ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas                         &
     2517     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas           &
     2518     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas                            &
     2519     &         ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas                              &
     2520     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas               &
     2521     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc)
     2522 
     2523       implicit none
     2524 
     2525#include "dimensions.h"
     2526
     2527!-------------------------------------------------------------------------
     2528! Vertical interpolation of TOGA-COARE forcing data onto mod_casel levels
     2529!-------------------------------------------------------------------------
     2530 
     2531       integer nlevmax
     2532       parameter (nlevmax=41)
     2533       integer nlev_cas,mxcalc
     2534!       real play(llm), plev_prof(nlevmax)
     2535!       real t_prof(nlevmax),q_prof(nlevmax)
     2536!       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
     2537!       real ht_prof(nlevmax),vt_prof(nlevmax)
     2538!       real hq_prof(nlevmax),vq_prof(nlevmax)
     2539 
     2540       real play(llm), plev_prof_cas(nlev_cas)
     2541       real t_prof_cas(nlev_cas),q_prof_cas(nlev_cas)
     2542       real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
     2543       real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas)
     2544       real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
     2545       real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
     2546       real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas),dtrad_prof_cas(nlev_cas)
     2547       real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
     2548 
     2549       real t_mod_cas(llm),q_mod_cas(llm)
     2550       real u_mod_cas(llm),v_mod_cas(llm)
     2551       real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm)
     2552       real du_mod_cas(llm),hu_mod_cas(llm),vu_mod_cas(llm)
     2553       real dv_mod_cas(llm),hv_mod_cas(llm),vv_mod_cas(llm)
     2554       real dt_mod_cas(llm),ht_mod_cas(llm),vt_mod_cas(llm),dtrad_mod_cas(llm)
     2555       real dq_mod_cas(llm),hq_mod_cas(llm),vq_mod_cas(llm)
     2556 
     2557       integer l,k,k1,k2
     2558       real frac,frac1,frac2,fact
     2559 
     2560       do l = 1, llm
     2561
     2562        if (play(l).ge.plev_prof_cas(nlev_cas)) then
     2563 
     2564        mxcalc=l
     2565         k1=0
     2566         k2=0
     2567
     2568         if (play(l).le.plev_prof_cas(1)) then
     2569
     2570         do k = 1, nlev_cas-1
     2571          if (play(l).le.plev_prof_cas(k).and. play(l).gt.plev_prof_cas(k+1)) then
     2572            k1=k
     2573            k2=k+1
     2574          endif
     2575         enddo
     2576
     2577         if (k1.eq.0 .or. k2.eq.0) then
     2578          write(*,*) 'PB! k1, k2 = ',k1,k2
     2579          write(*,*) 'l,play(l) = ',l,play(l)/100
     2580         do k = 1, nlev_cas-1
     2581          write(*,*) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/100
     2582         enddo
     2583         endif
     2584
     2585         frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))
     2586         t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1))
     2587         q_mod_cas(l)= q_prof_cas(k2) - frac*(q_prof_cas(k2)-q_prof_cas(k1))
     2588         u_mod_cas(l)= u_prof_cas(k2) - frac*(u_prof_cas(k2)-u_prof_cas(k1))
     2589         v_mod_cas(l)= v_prof_cas(k2) - frac*(v_prof_cas(k2)-v_prof_cas(k1))
     2590         ug_mod_cas(l)= ug_prof_cas(k2) - frac*(ug_prof_cas(k2)-ug_prof_cas(k1))
     2591         vg_mod_cas(l)= vg_prof_cas(k2) - frac*(vg_prof_cas(k2)-vg_prof_cas(k1))
     2592         w_mod_cas(l)= vitw_prof_cas(k2) - frac*(vitw_prof_cas(k2)-vitw_prof_cas(k1))
     2593         du_mod_cas(l)= du_prof_cas(k2) - frac*(du_prof_cas(k2)-du_prof_cas(k1))
     2594         hu_mod_cas(l)= hu_prof_cas(k2) - frac*(hu_prof_cas(k2)-hu_prof_cas(k1))
     2595         vu_mod_cas(l)= vu_prof_cas(k2) - frac*(vu_prof_cas(k2)-vu_prof_cas(k1))
     2596         dv_mod_cas(l)= dv_prof_cas(k2) - frac*(dv_prof_cas(k2)-dv_prof_cas(k1))
     2597         hv_mod_cas(l)= hv_prof_cas(k2) - frac*(hv_prof_cas(k2)-hv_prof_cas(k1))
     2598         vv_mod_cas(l)= vv_prof_cas(k2) - frac*(vv_prof_cas(k2)-vv_prof_cas(k1))
     2599         dt_mod_cas(l)= dt_prof_cas(k2) - frac*(dt_prof_cas(k2)-dt_prof_cas(k1))
     2600         ht_mod_cas(l)= ht_prof_cas(k2) - frac*(ht_prof_cas(k2)-ht_prof_cas(k1))
     2601         vt_mod_cas(l)= vt_prof_cas(k2) - frac*(vt_prof_cas(k2)-vt_prof_cas(k1))
     2602         dq_mod_cas(l)= dq_prof_cas(k2) - frac*(dq_prof_cas(k2)-dq_prof_cas(k1))
     2603         hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))
     2604         vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))
     2605     
     2606         else !play>plev_prof_cas(1)
     2607
     2608         k1=1
     2609         k2=2
     2610         frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))
     2611         frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))
     2612         t_mod_cas(l)= frac1*t_prof_cas(k1) - frac2*t_prof_cas(k2)
     2613         q_mod_cas(l)= frac1*q_prof_cas(k1) - frac2*q_prof_cas(k2)
     2614         u_mod_cas(l)= frac1*u_prof_cas(k1) - frac2*u_prof_cas(k2)
     2615         v_mod_cas(l)= frac1*v_prof_cas(k1) - frac2*v_prof_cas(k2)
     2616         ug_mod_cas(l)= frac1*ug_prof_cas(k1) - frac2*ug_prof_cas(k2)
     2617         vg_mod_cas(l)= frac1*vg_prof_cas(k1) - frac2*vg_prof_cas(k2)
     2618         w_mod_cas(l)= frac1*vitw_prof_cas(k1) - frac2*vitw_prof_cas(k2)
     2619         du_mod_cas(l)= frac1*du_prof_cas(k1) - frac2*du_prof_cas(k2)
     2620         hu_mod_cas(l)= frac1*hu_prof_cas(k1) - frac2*hu_prof_cas(k2)
     2621         vu_mod_cas(l)= frac1*vu_prof_cas(k1) - frac2*vu_prof_cas(k2)
     2622         dv_mod_cas(l)= frac1*dv_prof_cas(k1) - frac2*dv_prof_cas(k2)
     2623         hv_mod_cas(l)= frac1*hv_prof_cas(k1) - frac2*hv_prof_cas(k2)
     2624         vv_mod_cas(l)= frac1*vv_prof_cas(k1) - frac2*vv_prof_cas(k2)
     2625         dt_mod_cas(l)= frac1*dt_prof_cas(k1) - frac2*dt_prof_cas(k2)
     2626         ht_mod_cas(l)= frac1*ht_prof_cas(k1) - frac2*ht_prof_cas(k2)
     2627         vt_mod_cas(l)= frac1*vt_prof_cas(k1) - frac2*vt_prof_cas(k2)
     2628         dq_mod_cas(l)= frac1*dq_prof_cas(k1) - frac2*dq_prof_cas(k2)
     2629         hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2)
     2630         vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2)
     2631
     2632         endif ! play.le.plev_prof_cas(1)
     2633
     2634        else ! above max altitude of forcing file
     2635 
     2636!jyg
     2637         fact=20.*(plev_prof_cas(nlev_cas)-play(l))/plev_prof_cas(nlev_cas) !jyg
     2638         fact = max(fact,0.)                                           !jyg
     2639         fact = exp(-fact)                                             !jyg
     2640         t_mod_cas(l)= t_prof_cas(nlev_cas)                                   !jyg
     2641         q_mod_cas(l)= q_prof_cas(nlev_cas)*fact                              !jyg
     2642         u_mod_cas(l)= u_prof_cas(nlev_cas)*fact                              !jyg
     2643         v_mod_cas(l)= v_prof_cas(nlev_cas)*fact                              !jyg
     2644         ug_mod_cas(l)= ug_prof_cas(nlev_cas)*fact                              !jyg
     2645         vg_mod_cas(l)= vg_prof_cas(nlev_cas)*fact                              !jyg
     2646         w_mod_cas(l)= 0.0                                                 !jyg
     2647         du_mod_cas(l)= du_prof_cas(nlev_cas)*fact
     2648         hu_mod_cas(l)= hu_prof_cas(nlev_cas)*fact                            !jyg
     2649         vu_mod_cas(l)= vu_prof_cas(nlev_cas)*fact                            !jyg
     2650         dv_mod_cas(l)= dv_prof_cas(nlev_cas)*fact
     2651         hv_mod_cas(l)= hv_prof_cas(nlev_cas)*fact                            !jyg
     2652         vv_mod_cas(l)= vv_prof_cas(nlev_cas)*fact                            !jyg
     2653         dt_mod_cas(l)= dt_prof_cas(nlev_cas)
     2654         ht_mod_cas(l)= ht_prof_cas(nlev_cas)                                 !jyg
     2655         vt_mod_cas(l)= vt_prof_cas(nlev_cas)                                 !jyg
     2656         dq_mod_cas(l)= dq_prof_cas(nlev_cas)*fact
     2657         hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact                            !jyg
     2658         vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact                            !jyg
     2659 
     2660        endif ! play
     2661 
     2662       enddo ! l
     2663
     2664!       do l = 1,llm
     2665!       print *,'t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l) ',
     2666!     $        l,t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l)
     2667!       enddo
     2668 
     2669          return
     2670          end
     2671!*****************************************************************************
     2672!=====================================================================
    24332673       SUBROUTINE interp_dice_vertical(play,nlev_dice,nt_dice,plev_prof   &
    24342674     &         ,th_prof,qv_prof,u_prof,v_prof,o3_prof                     &

  • LMDZ5/branches/testing/libf/phylmd/1D_decl_cases.h

    r2160 r2220  
    240240        real u_profa(nlev_astex),v_profa(nlev_astex),w_profa(nlev_astex)
    241241        real tke_profa(nlev_astex)
    242 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    243 
     242
     243!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     244!Declarations specifiques au cas standard
     245
     246        real w_mod_cas(llm), t_mod_cas(llm),q_mod_cas(llm)
     247        real ug_mod_cas(llm),vg_mod_cas(llm)
     248        real u_mod_cas(llm),v_mod_cas(llm)
     249        real ht_mod_cas(llm),vt_mod_cas(llm),dt_mod_cas(llm),dtrad_mod_cas(llm)
     250        real hq_mod_cas(llm),vq_mod_cas(llm),dq_mod_cas(llm)
     251        real hu_mod_cas(llm),vu_mod_cas(llm),du_mod_cas(llm)
     252        real hv_mod_cas(llm),vv_mod_cas(llm),dv_mod_cas(llm)
     253!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     254

  • LMDZ5/branches/testing/libf/phylmd/1D_interp_cases.h

    r2160 r2220  
    597597      enddo
    598598      endif ! forcing_astex
    599 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    600 
     599
     600!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     601!---------------------------------------------------------------------
     602! Interpolation forcing standard case
     603!---------------------------------------------------------------------
     604      if (forcing_case) then
     605
     606        print*,                                                             &
     607     & '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/pdt_cas=',     &
     608     &    daytime,day1,(daytime-day1)*86400.,                               &
     609     &    (daytime-day1)*86400/pdt_cas
     610
     611! time interpolation:
     612        CALL interp_case_time(daytime,day1,annee_ref                        &
     613     &       ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas       &
     614     &       ,ts_cas,plev_cas,t_cas,q_cas,u_cas,v_cas,ug_cas,vg_cas    &
     615     &       ,vitw_cas,du_cas,hu_cas,vu_cas           &
     616     &       ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas             &
     617     &       ,dq_cas,hq_cas,vq_cas,lat_cas,sens_cas                 &
     618     &       ,ts_prof_cas,plev_prof_cas,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas         &
     619     &       ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas     &
     620     &       ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas       &
     621     &       ,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas,lat_prof_cas               &
     622     &       ,sens_prof_cas)
     623
     624             ts_cur = ts_prof_cas
     625             psurf=plev_prof_cas(1)
     626
     627! vertical interpolation:
     628      CALL interp_case_vertical(play,nlev_cas,plev_prof_cas            &
     629     &         ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas                         &
     630     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas           &
     631     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas           &
     632     &         ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas                              &
     633     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas               &
     634     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc)
     635
     636
     637!calcul de l'advection verticale a partir du omega
     638!Calcul des gradients verticaux
     639!initialisation
     640      d_t_z(:)=0.
     641      d_q_z(:)=0.
     642      d_u_z(:)=0.
     643      d_v_z(:)=0.
     644      d_t_dyn_z(:)=0.
     645      d_q_dyn_z(:)=0.
     646      d_u_dyn_z(:)=0.
     647      d_v_dyn_z(:)=0.
     648      DO l=2,llm-1
     649       d_t_z(l)=(temp(l+1)-temp(l-1))/(play(l+1)-play(l-1))
     650       d_q_z(l)=(q(l+1,1)-q(l-1,1))/(play(l+1)-play(l-1))
     651       d_u_z(l)=(u(l+1)-u(l-1))/(play(l+1)-play(l-1))
     652       d_v_z(l)=(v(l+1)-v(l-1))/(play(l+1)-play(l-1))
     653      ENDDO
     654      d_t_z(1)=d_t_z(2)
     655      d_q_z(1)=d_q_z(2)
     656      d_u_z(1)=d_u_z(2)
     657      d_v_z(1)=d_v_z(2)
     658      d_t_z(llm)=d_t_z(llm-1)
     659      d_q_z(llm)=d_q_z(llm-1)
     660      d_u_z(llm)=d_u_z(llm-1)
     661      d_v_z(llm)=d_v_z(llm-1)
     662
     663!Calcul de l advection verticale
     664      d_t_dyn_z(:)=w_mod_cas(:)*d_t_z(:)
     665      d_q_dyn_z(:)=w_mod_cas(:)*d_q_z(:)
     666      d_u_dyn_z(:)=w_mod_cas(:)*d_u_z(:)
     667      d_v_dyn_z(:)=w_mod_cas(:)*d_v_z(:)
     668
     669!wind nudging
     670      if (nudge_u.gt.0.) then
     671        do l=1,llm
     672           u(l)=u(l)+timestep*(u_mod_cas(l)-u(l))/(nudge_u)
     673        enddo
     674      else
     675        do l=1,llm
     676        u(l) = u_mod_cas(l)
     677        enddo
     678      endif
     679
     680      if (nudge_v.gt.0.) then
     681        do l=1,llm
     682           v(l)=v(l)+timestep*(v_mod_cas(l)-v(l))/(nudge_v)
     683        enddo
     684      else
     685        do l=1,llm
     686        v(l) = v_mod_cas(l)
     687        enddo
     688      endif
     689
     690      if (nudge_w.gt.0.) then
     691        do l=1,llm
     692           w(l)=w(l)+timestep*(w_mod_cas(l)-w(l))/(nudge_w)
     693        enddo
     694      else
     695        do l=1,llm
     696        w(l) = w_mod_cas(l)
     697        enddo
     698      endif
     699
     700!nudging of q and temp
     701      if (nudge_t.gt.0.) then
     702        do l=1,llm
     703           temp(l)=temp(l)+timestep*(t_mod_cas(l)-temp(l))/(nudge_t)
     704        enddo
     705      endif
     706      if (nudge_q.gt.0.) then
     707        do l=1,llm
     708           q(l,1)=q(l,1)+timestep*(q_mod_cas(l)-q(l,1))/(nudge_q)
     709        enddo
     710      endif
     711
     712      do l = 1, llm
     713       omega(l) = w_mod_cas(l)
     714       omega2(l)= omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
     715       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
     716
     717!calcul advection
     718        if ((tend_u.eq.1).and.(tend_w.eq.0)) then
     719           d_u_adv(l)=du_mod_cas(l)
     720        else if ((tend_u.eq.1).and.(tend_w.eq.1)) then
     721           d_u_adv(l)=hu_mod_cas(l)-d_u_dyn_z(l)
     722        endif
     723
     724        if ((tend_v.eq.1).and.(tend_w.eq.0)) then
     725           d_v_adv(l)=dv_mod_cas(l)
     726        else if ((tend_v.eq.1).and.(tend_w.eq.1)) then
     727           d_v_adv(l)=hv_mod_cas(l)-d_v_dyn_z(l)
     728        endif
     729
     730        if ((tend_t.eq.1).and.(tend_w.eq.0)) then
     731!           d_th_adv(l)=alpha*omega(l)/rcpd+dt_mod_cas(l)
     732           d_th_adv(l)=alpha*omega(l)/rcpd-dt_mod_cas(l)
     733        else if ((tend_t.eq.1).and.(tend_w.eq.1)) then
     734!           d_th_adv(l)=alpha*omega(l)/rcpd+ht_mod_cas(l)-d_t_dyn_z(l)
     735           d_th_adv(l)=alpha*omega(l)/rcpd-ht_mod_cas(l)-d_t_dyn_z(l)
     736        endif
     737
     738        if ((tend_q.eq.1).and.(tend_w.eq.0)) then
     739!           d_q_adv(l,1)=dq_mod_cas(l)
     740           d_q_adv(l,1)=-1*dq_mod_cas(l)
     741        else if ((tend_q.eq.1).and.(tend_w.eq.1)) then
     742!           d_q_adv(l,1)=hq_mod_cas(l)-d_q_dyn_z(l)
     743           d_q_adv(l,1)=-1*hq_mod_cas(l)-d_q_dyn_z(l)
     744        endif
     745         
     746        if (tend_rayo.eq.1) then
     747           dt_cooling(l) = dtrad_mod_cas(l)
     748        else
     749           dt_cooling(l) = 0.0
     750        endif
     751      enddo
     752
     753      endif ! forcing_case
     754
     755
     756!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     757

  • LMDZ5/branches/testing/libf/phylmd/1D_read_forc_cases.h

    r2160 r2220  
    720720
    721721      endif ! forcing_astex
    722 
     722!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     723!---------------------------------------------------------------------
     724! Forcing from standard case :
     725!---------------------------------------------------------------------
     726
     727      if (forcing_case) then
     728
     729         write(*,*),'avant call read_1D_cas'
     730         call read_1D_cas
     731         write(*,*) 'Forcing read'
     732
     733!Time interpolation for initial conditions using TOGA interpolation routine
     734         write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1   
     735      CALL interp_case_time(day,day1,annee_ref                &
     736     &         ,year_ini_cas,day_ju_ini_cas,nt_cas,pdt_cas,nlev_cas       &
     737     &         ,ts_cas,plev_cas,t_cas,q_cas,u_cas,v_cas               &
     738     &         ,ug_cas,vg_cas,vitw_cas,du_cas,hu_cas,vu_cas           &
     739     &         ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas             &
     740     &         ,dq_cas,hq_cas,vq_cas,lat_cas,sens_cas                 &
     741     &         ,ts_prof_cas,plev_prof_cas,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas         &
     742     &         ,ug_prof_cas,vg_prof_cas,vitw_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas     &
     743     &         ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas       &
     744     &         ,dq_prof_cas,hq_prof_cas,vq_prof_cas,lat_prof_cas,sens_prof_cas)
     745
     746! vertical interpolation using TOGA interpolation routine:
     747!      write(*,*)'avant interp vert', t_prof
     748      CALL interp_case_vertical(play,nlev_cas,plev_prof_cas            &
     749     &         ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas    &
     750     &         ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas           &
     751     &         ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas           &
     752     &         ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas           &
     753     &         ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas               &
     754     &         ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc)
     755!       write(*,*) 'Profil initial forcing case interpole',t_mod
     756
     757! initial and boundary conditions :
     758!      tsurf = ts_prof_cas
     759      ts_cur = ts_prof_cas
     760      psurf=plev_prof_cas(1)
     761      write(*,*) 'SST initiale: ',tsurf
     762      do l = 1, llm
     763       temp(l) = t_mod_cas(l)
     764       q(l,1) = q_mod_cas(l)
     765       q(l,2) = 0.0
     766       u(l) = u_mod_cas(l)
     767       v(l) = v_mod_cas(l)
     768       omega(l) = w_mod_cas(l)
     769       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
     770
     771       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
     772!on applique le forcage total au premier pas de temps
     773!attention: signe different de toga
     774       d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod_cas(l)+vt_mod_cas(l))
     775       d_q_adv(l,1) = (hq_mod_cas(l)+vq_mod_cas(l))
     776       d_q_adv(l,2) = 0.0
     777       d_u_adv(l) = (hu_mod_cas(l)+vu_mod_cas(l))
     778       d_u_adv(l) = (hv_mod_cas(l)+vv_mod_cas(l))
     779      enddo     
     780       
     781      endif !forcing_case
     782!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  • LMDZ5/branches/testing/libf/phylmd/YOMCST2.h

    r1910 r2220  
    33      REAL  gammas, alphas, betas, Fmax, qqa1, qqa2, qqa3, scut
    44      REAL  Qcoef1max,Qcoef2max,Supcrit1,Supcrit2
     5      REAL coef_clos_ls
    56!
    67      COMMON/YOMCST2/gammas,    alphas, betas, Fmax, scut,              &
    … …  
    89     &               Qcoef1max,Qcoef2max,                               &
    910     &               Supcrit1, Supcrit2,                                &
    10      &               choice,iflag_mix
     11     &               choice,iflag_mix,coef_clos_ls
    1112!$OMP THREADPRIVATE(/YOMCST2/)
    1213!    --------------------------------------------------------------------

  • LMDZ5/branches/testing/libf/phylmd/calcratqs.F90

    r1910 r2220  
    11SUBROUTINE calcratqs(klon,klev,prt_level,lunout,       &
    2            iflag_ratqs,iflag_con,iflag_cldcon,pdtphys, &
     2           iflag_ratqs,iflag_con,iflag_cldth,pdtphys, &
    33           ratqsbas,ratqshaut,tau_ratqs,fact_cldcon,   &
    44           ptconv,ptconvth,clwcon0th, rnebcon0th,      &
    … …  
    1919! Input
    2020integer,intent(in) :: klon,klev,prt_level,lunout
    21 integer,intent(in) :: iflag_con,iflag_cldcon,iflag_ratqs
     21integer,intent(in) :: iflag_con,iflag_cldth,iflag_ratqs
    2222real,intent(in) :: pdtphys,ratqsbas,ratqshaut,fact_cldcon,tau_ratqs
    2323real, dimension(klon,klev+1),intent(in) :: paprs
    … …  
    4343!   ----------------
    4444!   on ecrase le tableau ratqsc calcule par clouds_gno
    45       if (iflag_cldcon.eq.1) then
     45      if (iflag_cldth.eq.1) then
    4646         do k=1,klev
    4747         do i=1,klon
    … …  
    5858!  par nversion de la fonction log normale
    5959!-----------------------------------------------------------------------
    60       else if (iflag_cldcon.eq.4) then
     60      else if (iflag_cldth.eq.4) then
    6161         ptconvth(:,:)=.false.
    6262         ratqsc(:,:)=0.
    … …  
    136136!  -----------
    137137
    138       if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2.or.iflag_cldcon.eq.4) then
     138      if (iflag_cldth.eq.1 .or.iflag_cldth.eq.2.or.iflag_cldth.eq.4) then
    139139
    140140! On ajoute une constante au ratqsc*2 pour tenir compte de
    … …  
    165165!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    166166         ratqs(:,:)=max(ratqs(:,:),ratqss(:,:))
    167       else if (iflag_cldcon<=6) then
     167      else if (iflag_cldth<=6) then
    168168!   on ne prend que le ratqs stable pour fisrtilp
    169169         ratqs(:,:)=ratqss(:,:)
    … …  
    174174             do i=1,klon
    175175                if (ratqsc(i,k).gt.1.e-10) then
    176                    ratqs(i,k)=ratqs(i,k)*zfratqs2+(iflag_cldcon/100.)*ratqsc(i,k)*(1.-zfratqs2)
     176                   ratqs(i,k)=ratqs(i,k)*zfratqs2+(iflag_cldth/100.)*ratqsc(i,k)*(1.-zfratqs2)
    177177                endif
    178178                ratqs(i,k)=min(ratqs(i,k)*zfratqs1+ratqss(i,k)*(1.-zfratqs1),0.5)

  • LMDZ5/branches/testing/libf/phylmd/change_srf_frac_mod.F90

    r2187 r2220  
    2323
    2424    USE dimphy
    25     USE surface_data, ONLY : type_ocean
     25    USE surface_data, ONLY : type_ocean,version_ocean
    2626    USE limit_read_mod
    2727    USE pbl_surface_mod, ONLY : pbl_surface_newfrac
    2828    USE cpl_mod, ONLY : cpl_receive_frac
    29     USE ocean_slab_mod, ONLY : ocean_slab_frac
     29    USE ocean_slab_mod, ONLY : fsic, ocean_slab_frac
    3030    USE indice_sol_mod
    3131
    … …  
    146146          pctsrf(:,is_sic) = 0.
    147147       END WHERE
     148! Send fractions back to slab ocean if needed
     149       IF (type_ocean == 'slab'.AND. version_ocean.NE.'sicINT') THEN
     150           WHERE (1.-zmasq(:)>EPSFRA)
     151               fsic(:)=pctsrf(:,is_sic)/(1.-zmasq(:))
     152           END WHERE
     153       END IF
    148154
    149155!****************************************************************************************

  • LMDZ5/branches/testing/libf/phylmd/clift.F90

    r1999 r2220  
    33
    44SUBROUTINE clift(p, t, rr, rs, plcl, dplcldt, dplcldq)
     5IMPLICIT NONE
    56  ! ***************************************************************
    67  ! *                                                             *
    … …  
    4142
    4243  include "YOMCST.h"
     44  real :: p,t,rr,rs,plcl,dplcldt,dplcldq,cpd,cpv,cl,cpvmcl,eps,alv0,a,b
     45  real :: rh,chi,alv
    4346
    4447  cpd = rcpd

  • LMDZ5/branches/testing/libf/phylmd/compar1d.h

    r2187 r2220  
    33!
    44      integer :: forcing_type
     5      integer :: tend_u,tend_v,tend_w,tend_t,tend_q,tend_rayo
     6      real :: nudge_u,nudge_v,nudge_w,nudge_t,nudge_q
    57      integer :: iflag_nudge
    68      real :: nat_surf
    … …  
    3234     & qsol,qsurf,psurf,zsurf,albedo,time,time_ini,xlat,xlon,airefi,    &
    3335     & wtsurf,wqsurf,restart_runoff,xagesno,qsolinp,zpicinp,            &
    34      & forcing_type,                                                    &
     36     & forcing_type,tend_u,tend_v,tend_w,tend_t,tend_q,tend_rayo,       &
     37     & nudge_u,nudge_v,nudge_w,nudge_t,nudge_q,                         &
    3538     & iflag_nudge,                                                     &
    3639     & restart,ok_old_disvert

  • LMDZ5/branches/testing/libf/phylmd/concvl.F90

    r2056 r2220  
    1515                  dd_t, dd_q, lalim_conv, wght_th, &                 ! RomP
    1616                  evap, ep, epmlmMm, eplaMm, &                       ! RomP
    17                   wdtrainA, wdtrainM, wght)                          ! RomP+RL
     17                  wdtrainA, wdtrainM, wght, qtc, sigt, &
     18                  tau_cld_cv, coefw_cld_cv)                           ! RomP+RL, AJ
    1819!RomP <<<
    1920! **************************************************************
    … …  
    2930  USE dimphy
    3031  USE infotrac, ONLY: nbtr
     32  USE phys_local_var_mod, ONLY: omega
    3133  IMPLICIT NONE
    3234! ======================================================================
    … …  
    135137  REAL dplcldt(klon), dplcldr(klon)
    136138  REAL qcondc(klon, klev)
     139  REAL qtc(klon, klev)
     140  REAL sigt(klon, klev)
    137141  REAL wd(klon)
    138142  REAL plim1(klon), plim2(klon), asupmax(klon, klev)
    … …  
    141145  REAL sigd(klon)
    142146  REAL zx_t, zdelta, zx_qs, zcor
     147  REAL tau_cld_cv, coefw_cld_cv
    143148
    144149!   INTEGER iflag_mix
    … …  
    203208!$OMP THREADPRIVATE(itap, igout)
    204209
     210
    205211  include "YOMCST.h"
    206212  include "YOMCST2.h"
    … …  
    398404                    t, q, qs, t_wake, q_wake, qs_wake, s_wake, u, v, tra, &
    399405                    em_p, em_ph, &
    400                     Ale, Alp, &
     406                    Ale, Alp, omega, &
    401407                    em_sig1feed, em_sig2feed, em_wght, &
    402408                    iflag, d_t, d_q, d_u, d_v, d_tra, rain, kbas, ktop, &
    … …  
    411417                    da, phi, mp, phi2, d1a, dam, sij, wght, &           ! RomP+RL
    412418                    clw, elij, evap, ep, epmlmMm, eplaMm, &             ! RomP+RL
    413                     wdtrainA, wdtrainM)                                 ! RomP
     419                    wdtrainA, wdtrainM, qtc, sigt, &
     420                    tau_cld_cv, coefw_cld_cv)                           ! RomP,AJ
    414421!AC!+!RomP+jyg
    415422  END IF

  • LMDZ5/branches/testing/libf/phylmd/conf_phys_m.F90

    r2187 r2220  
    1515       solarlong0,seuil_inversion, &
    1616       fact_cldcon, facttemps,ok_newmicro,iflag_radia,&
    17        iflag_cldcon, &
     17       iflag_cldth, &
    1818       iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
    1919       ok_ade, ok_aie, ok_cdnc, aerosol_couple, &
    … …  
    8181    REAL                 :: bl95_b0, bl95_b1
    8282    real                 :: fact_cldcon, facttemps,ratqsbas,ratqshaut,tau_ratqs
    83     integer              :: iflag_cldcon
     83    integer              :: iflag_cldth
    8484    integer              :: iflag_ratqs
    8585
    … …  
    101101    REAL,SAVE           :: freq_COSP_omp
    102102    real,SAVE           :: fact_cldcon_omp, facttemps_omp,ratqsbas_omp
     103    real,SAVE           :: tau_cld_cv_omp, coefw_cld_cv_omp
     104    integer,SAVE        :: iflag_cld_cv_omp
     105
     106
    103107    real,SAVE           :: ratqshaut_omp
    104108    real,SAVE           :: tau_ratqs_omp
    … …  
    107111    integer,SAVE        :: iflag_rrtm_omp
    108112    integer,SAVE        :: NSW_omp
    109     integer,SAVE        :: iflag_cldcon_omp, ip_ebil_phy_omp
     113    integer,SAVE        :: iflag_cldth_omp, ip_ebil_phy_omp
    110114    integer,SAVE        :: iflag_ratqs_omp
    111115
    … …  
    132136    integer,SAVE :: iflag_coupl_omp,iflag_clos_omp,iflag_wake_omp
    133137    integer,SAVE :: iflag_cvl_sigd_omp
     138    REAL, SAVE :: coef_clos_ls_omp
    134139    REAL, SAVE :: supcrit1_omp, supcrit2_omp
    135140    INTEGER, SAVE :: iflag_mix_omp
    … …  
    886891
    887892    !
    888     !Config Key  = iflag_cldcon
     893    !Config Key  = iflag_cldth
    889894    !Config Desc = 
    890895    !Config Def  = 1
    891896    !Config Help =
    892897    !
    893     iflag_cldcon_omp = 1
    894     call getin('iflag_cldcon',iflag_cldcon_omp)
     898    iflag_cldth_omp = 1
     899! On lit deux fois avec l'ancien et le nouveau nom
     900! pour assurer une retrocompatiblite.
     901! A abandonner un jour
     902    call getin('iflag_cldcon',iflag_cldth_omp)
     903    call getin('iflag_cldth',iflag_cldth_omp)
     904
     905    !
     906    !Config Key  = iflag_cld_cv
     907    !Config Desc =
     908    !Config Def  = 1
     909    !Config Help =
     910    !
     911    iflag_cld_cv_omp = 1
     912    call getin('iflag_cld_cv',iflag_cld_cv_omp)
     913
     914    !
     915    !Config Key  = tau_cld_cv
     916    !Config Desc =
     917    !Config Def  = 10.
     918    !Config Help =
     919    !
     920    tau_cld_cv_omp = 10.
     921    call getin('tau_cld_cv',tau_cld_cv_omp)
     922
     923    !
     924    !Config Key  = coefw_cld_cv
     925    !Config Desc =
     926    !Config Def  = 0.1
     927    !Config Help =
     928    !
     929    coefw_cld_cv_omp = 0.1
     930    call getin('coefw_cld_cv',coefw_cld_cv_omp)
     931
     932
     933
    895934
    896935    !
    … …  
    13431382    iflag_clos_omp = 1
    13441383    call getin('iflag_clos',iflag_clos_omp)
     1384    !
     1385    !Config Key  = coef_clos_ls
     1386    !Config Desc = 
     1387    !Config Def  = 0
     1388    !Config Help =
     1389    !
     1390    coef_clos_ls_omp = 0.
     1391    call getin('coef_clos_ls',coef_clos_ls_omp)
     1392
    13451393    !
    13461394    !Config Key  = iflag_cvl_sigd
    … …  
    19251973    iflag_rrtm = iflag_rrtm_omp
    19261974    NSW = NSW_omp
    1927     iflag_cldcon = iflag_cldcon_omp
     1975    iflag_cldth = iflag_cldth_omp
     1976    iflag_cld_cv = iflag_cld_cv_omp
     1977    tau_cld_cv = tau_cld_cv_omp
     1978    coefw_cld_cv = coefw_cld_cv_omp
    19281979    iflag_ratqs = iflag_ratqs_omp
    19291980    ip_ebil_phy = ip_ebil_phy_omp
    … …  
    19461997    iflag_clos = iflag_clos_omp
    19471998    iflag_wake = iflag_wake_omp
     1999    coef_clos_ls = coef_clos_ls_omp
    19482000    alp_offset = alp_offset_omp
    19492001    iflag_cvl_sigd = iflag_cvl_sigd_omp
    … …  
    20762128    write(lunout,*)' reevap_ice = ', reevap_ice
    20772129    write(lunout,*)' iflag_pdf = ', iflag_pdf
    2078     write(lunout,*)' iflag_cldcon = ', iflag_cldcon
     2130    write(lunout,*)' iflag_cldth = ', iflag_cldth
     2131    write(lunout,*)' iflag_cld_cv = ', iflag_cld_cv
     2132    write(lunout,*)' tau_cld_cv = ', tau_cld_cv
     2133    write(lunout,*)' coefw_cld_cv = ', coefw_cld_cv
    20792134    write(lunout,*)' iflag_radia = ', iflag_radia
    20802135    write(lunout,*)' iflag_rrtm = ', iflag_rrtm
    … …  
    21352190    write(lunout,*)' iflag_thermals_closure = ', iflag_thermals_closure
    21362191    write(lunout,*)' iflag_clos = ', iflag_clos
     2192    write(lunout,*)' coef_clos_ls = ', coef_clos_ls
    21372193    write(lunout,*)' type_run = ',type_run
    21382194    write(lunout,*)' ok_cosp = ',ok_cosp

  • LMDZ5/branches/testing/libf/phylmd/convect3.F90

    r1999 r2220  
    1717  USE dimphy
    1818  USE infotrac, ONLY: nbtr
    19 
     19  IMPLICIT NONE
    2020  include "dimensions.h"
    2121  INTEGER na
    … …  
    7373
    7474
    75 
     75  REAL :: cpv,cl,cpvmcl,eps,alv0,rdcp,pbcrit,ptcrit,sigd,spfac
     76  REAL :: tau,beta,alpha,dtcrit,dtovsh,ahm,rm,um,vm,dphinv
     77  REAL :: a2,x,tvx,tvy,plcl,pden,dpbase,tvpbase,tvbase,tdif
     78  REAL :: ath1,ath,delti,deltap,dcape,dlnp,sigold,dtmin,fac,w
     79  REAL :: amu,rti,cpd,bf2,anum,denom,dei,altem,cwat,stemp,qp
     80  REAL :: scrit,alt,smax,asij,wgh,sjmax,sjmin,smid,delp,delm
     81  REAL :: asum,bsum,csum,wflux,tinv,wdtrain,awat,afac,afac1,afac2
     82  REAL :: bfac,pr1,pr2,sigt,b6,c6,revap,tevap,delth,amfac,amp2
     83  REAL :: xf,tf,af,bf,fac2,ur,sru,d,ampmax,dpinv,am,amde,cpinv
     84  REAL :: amp1,ad,rat,ax,bx,cx,dx,ex,dsum
     85  INTEGER :: nk,i,j,nopt,jn,k,im,jm,n
    7686
    7787  REAL dnwd0(nd) !  precipitation driven unsaturated downdraft flux

  • LMDZ5/branches/testing/libf/phylmd/cv3_inicp.F90

    r1999 r2220  
    99  ! modified by :                                               *
    1010  ! **************************************************************
    11 
     11  IMPLICIT NONE
    1212  include "YOMCST2.h"
    1313
    … …  
    1919
    2020  REAL qcoef1, qcoef2, qff, qfff, qmix, rmix, qmix1, rmix1, qmix2, rmix2, f
     21  REAL :: sumcoef,sigma,aire,pdf,mu,df,ff
    2122
    2223  qcoef1(f) = tanh(f/gammas)

  • LMDZ5/branches/testing/libf/phylmd/cv3_inip.F90

    r1999 r2220  
    11SUBROUTINE cv3_inip()
    2   ! **************************************************************
    3   ! *
    4   ! CV3_INIP Lecture des choix de lois de probabilité de mélange*
    5   ! et calcul de leurs coefficients normalisés.        *
    6   ! *
    7   ! written by   : Jean-Yves Grandpeix, 06/06/2006, 19.39.27    *
    8   ! modified by :                                               *
    9   ! **************************************************************
     2  ! *******************************************************************
     3  ! *                                                                  *
     4  ! CV3_INIP Input = choice of mixing probability laws                 *
     5  !          Output = normalized coefficients of the probability laws. *
     6  ! *                                                                  *
     7  ! written by   : Jean-Yves Grandpeix, 06/06/2006, 19.39.27           *
     8  ! modified by :                                                      *
     9  ! *******************************************************************
     10!
     11!----------------------------------------------
     12!  INPUT (from Common YOMCST2 in "YOMCST2.h") :
     13! iflag_mix
     14! gammas   
     15! alphas
     16! betas
     17! Fmax
     18! scut
     19!
     20!----------------------------------------------
     21!  INPUT/OUTPUT (from and to Common YOMCST2 in "YOMCST2.h") :
     22! qqa1
     23! qqa2
     24!
     25!----------------------------------------------
     26!  OUTPUT (to Common YOMCST2 in "YOMCST2.h") :
     27! Qcoef1max
     28! Qcoef2max
     29!
     30!----------------------------------------------
     31
     32  IMPLICIT NONE
    1033
    1134  include "YOMCST2.h"
    1235
    13   ! INTEGER iflag_mix
    1436  include 'iniprint.h'
    1537
    16   CHARACTER (LEN=20) :: modname = 'cv3_inip'
    17   CHARACTER (LEN=80) :: abort_message
     38!----------------------------------------------
     39! Local variables :
     40   CHARACTER (LEN=20) :: modname = 'cv3_inip'
     41   CHARACTER (LEN=80) :: abort_message
     42
     43   REAL               :: sumcoef
     44   REAL               :: sigma, aire, pdf, mu, df
     45   REAL               :: ff
    1846
    1947

  • LMDZ5/branches/testing/libf/phylmd/cv3_routines.F90

    r2160 r2220  
    11
    22! $Id$
     3
    34
    45
    … …  
    27632764                     cbmf, upwd, dnwd, dnwd0, ma, mip, &
    27642765                     tls, tps, qcondc, wd, &
    2765                      ftd, fqd)
     2766                     ftd, fqd, qnk, qtc, sigt, tau_cld_cv, coefw_cld_cv)
    27662767
    27672768  IMPLICIT NONE
    … …  
    28002801!input/output:
    28012802      INTEGER iflag(nloc)
     2803      REAL,INTENT(in) :: tau_cld_cv, coefw_cld_cv
    28022804!
    28032805!outputs:
    … …  
    28112813      REAL tls(nloc, nd), tps(nloc, nd)
    28122814      REAL qcondc(nloc, nd) ! cld
     2815      REAL qtc(nloc,nd), sigt(nloc,nd) ! cld
    28132816      REAL wd(nloc) ! gust
    28142817      REAL cbmf(nloc)
    … …  
    28302833      REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld
    28312834      REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld
    2832 
     2835      REAL sument(nloc), sigment(nloc,nd), qtment(nloc,nd) ! cld
     2836      REAL qnk(nloc)
    28332837      REAL sumdq !jyg
    28342838!
    … …  
    28612865      qcondc(il, i) = 0.0 ! cld
    28622866      qcond(il, i) = 0.0 ! cld
     2867      qtc(il, i) = 0.0 ! cld
     2868      qtment(il, i) = 0.0 ! cld
     2869      sigment(il, i) = 0.0 ! cld
     2870      sigt(il, i) = 0.0 ! cld
    28632871      nqcond(il, i) = 0.0 ! cld
    28642872    END DO
    … …  
    32343242! (saturated updrafts resulting from mixing)                                   ! cld
    32353243          qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat(il))                ! cld
    3236           nqcond(il, i) = nqcond(il, i) + 1. ! cld
     3244          qtment(il, i) = qtment(il, i) + qent(il,k,i)                         ! cld
     3245          nqcond(il, i) = nqcond(il, i) + 1.                                   ! cld
    32373246        END IF ! i
    32383247      END DO
    … …  
    33103319! (saturated downdrafts resulting from mixing)                                 ! cld
    33113320          qcond(il, i) = qcond(il, i) + elij(il, k, i)                         ! cld
    3312           nqcond(il, i) = nqcond(il, i) + 1. ! cld
     3321          qtment(il, i) = qent(il,k,i) + qtment(il,i)                          ! cld
     3322          nqcond(il, i) = nqcond(il, i) + 1.                                   ! cld
    33133323        END IF ! cld
    33143324      END DO ! cld
    … …  
    33193329      IF (i<=inb(il) .AND. nent(il,i)==0 .AND. iflag(il)<=1) THEN              ! cld
    33203330        qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i)                 ! cld
     3331        qtment(il, i) = qent(il,k,i) + qtment(il,i)                          ! cld
    33213332        nqcond(il, i) = nqcond(il, i) + 1.                                     ! cld
    33223333      END IF                                                                   ! cld
    … …  
    33263337      IF (i<=inb(il) .AND. nqcond(il,i)/=0 .AND. iflag(il)<=1) THEN            ! cld
    33273338        qcond(il, i) = qcond(il, i)/nqcond(il, i)                              ! cld
     3339        qtment(il, i) = qtment(il,i)/nqcond(il, i)                             ! cld
    33283340      END IF                                                                   ! cld
    33293341    END DO
    … …  
    37883800      END IF                                                         ! cld
    37893801    END DO                                                           ! cld
    3790   END DO                                                             ! cld
    3791 
    3792   DO i = 1, nl                                                       ! cld
     3802  END DO 
     3803                                                           ! cld
     3804  DO i = 1, nl 
     3805
     3806! 14/01/15 AJ je remets les parties manquantes cf JYG
     3807! Initialize sument to 0
     3808
     3809    DO il = 1,ncum
     3810     sument(il) = 0.
     3811    ENDDO
     3812
     3813! Sum mixed mass fluxes in sument
     3814
     3815    DO k = 1,nl
     3816      DO il = 1,ncum
     3817        IF  (k<=inb(il) .AND. i<=inb(il) .AND. iflag(il)<=1) THEN   ! cld
     3818          sument(il) =sument(il) + abs(ment(il,k,i))
     3819        ENDIF
     3820      ENDDO     ! il
     3821    ENDDO       ! k
     3822
     3823! 14/01/15 AJ delta n'a rien à faire là...                                                 
    37933824    DO il = 1, ncum                                                  ! cld
    37943825      IF (wa(il,i)>0.0 .AND. iflag(il)<=1) &                         ! cld
    3795         siga(il, i) = mac(il, i)/wa(il, i) &                         ! cld
    3796         *rrd*tvp(il, i)/p(il, i)/100./delta                          ! cld
     3826        siga(il, i) = mac(il, i)/(coefw_cld_cv*wa(il, i)) &          ! cld
     3827        *rrd*tvp(il, i)/p(il, i)/100.                                ! cld
     3828
    37973829      siga(il, i) = min(siga(il,i), 1.0)                             ! cld
    3798 ! IM cf. FH                                                         
     3830
     3831! IM cf. FH
     3832! 14/01/15 AJ ne correspond pas à ce qui a été codé par JYG et SB           
     3833                                                         
    37993834      IF (iflag_clw==0) THEN                                         ! cld
    38003835        qcondc(il, i) = siga(il, i)*clw(il, i)*(1.-ep(il,i)) &       ! cld
    38013836          +(1.-siga(il,i))*qcond(il, i)                              ! cld
     3837
     3838
     3839        sigment(il,i)=sument(il)*tau_cld_cv/(ph(il,i)-ph(il,i+1))    ! cld
     3840        sigment(il, i) = min(1.e-4+sigment(il,i), 1.0 - siga(il,i))  ! cld
     3841        qtc(il, i) = (siga(il,i)*qnk(il)+sigment(il,i)*qtment(il,i)) & ! cld
     3842                     /(siga(il,i)+sigment(il,i))                     ! cld
     3843        sigt(il,i) = sigment(il, i) + siga(il, i)
     3844
     3845!        qtc(il, i) = siga(il,i)*qnk(il)+(1.-siga(il,i))*qtment(il,i) ! cld
     3846!     print*,'BIGAUSSIAN CONV',siga(il,i),sigment(il,i),qtc(il,i) 
     3847               
    38023848      ELSE IF (iflag_clw==1) THEN                                    ! cld
    38033849        qcondc(il, i) = qcond(il, i)                                 ! cld
     3850        qtc(il,i) = qtment(il,i)                                     ! cld
    38043851      END IF                                                         ! cld
    38053852

  • LMDZ5/branches/testing/libf/phylmd/cv3a_compress.F90

    r1999 r2220  
    44    th1_wake, tra1, h1, lv1, lf1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, &
    55    h1_wake, lv1_wake, lf1_wake, cpn1_wake, tv1_wake, sig1, w01, ptop21, &
    6     ale1, alp1, iflag, nk, icb, icbs, plcl, tnk, qnk, gznk, hnk, unk, vnk, &
     6    ale1, alp1, omega1, iflag, nk, icb, icbs, plcl, tnk, qnk, gznk, hnk, unk, vnk, &
    77    wghti, pbase, buoybase, t, q, qs, t_wake, q_wake, qs_wake, s_wake, u, v, &
    88    gz, th, th_wake, tra, h, lv, lf, cpn, p, ph, tv, tp, tvp, clw, h_wake, &
    9     lv_wake, lf_wake, cpn_wake, tv_wake, sig, w0, ptop2, ale, alp)
     9    lv_wake, lf_wake, cpn_wake, tv_wake, sig, w0, ptop2, ale, alp, omega)
    1010  ! **************************************************************
    1111  ! *
    … …  
    4040  REAL sig1(len, nd), w01(len, nd), ptop21(len)
    4141  REAL ale1(len), alp1(len)
     42  REAL omega1(len,nd)
    4243
    4344  ! outputs:
    … …  
    6061  REAL sig(len, nd), w0(len, nd), ptop2(len)
    6162  REAL ale(len), alp(len)
     63  REAL omega(len,nd)
    6264
    6365  ! local variables:
    … …  
    102104        sig(nn, k) = sig1(i, k)
    103105        w0(nn, k) = w01(i, k)
     106        omega(nn, k) = omega1(i, k)
    104107      END IF
    105108    END DO

  • LMDZ5/branches/testing/libf/phylmd/cv3a_uncompress.F90

    r1999 r2220  
    66    , clw, elij, evap, ep, epmlmmm, eplamm & ! RomP
    77    , wdtraina, wdtrainm &         ! RomP
     8    , qtc, sigt          &
    89
    910    , iflag1, kbas1, ktop1, precip1, cbmf1, plcl1, plfc1, wbeff1, sig1, w01, &
    … …  
    1314    , da1, phi1, mp1, phi21, d1a1, dam1, sigij1 & ! RomP+AC+jyg
    1415    , clw1, elij1, evap1, ep1, epmlmmm1, eplamm1 & ! RomP
    15     , wdtraina1, wdtrainm1) ! RomP
     16    , wdtraina1, wdtrainm1 & ! RomP
     17    , qtc1, sigt1)
    1618
    1719  ! **************************************************************
    … …  
    5658  REAL evap(nloc, nd), ep(nloc, nd) !RomP
    5759  REAL epmlmmm(nloc, nd, nd), eplamm(nloc, nd) !RomP+jyg
     60  REAL qtc(nloc, nd), sigt(nloc, nd) !RomP
    5861  REAL wdtraina(nloc, nd), wdtrainm(nloc, nd) !RomP
    5962
    … …  
    8487  REAL evap1(len, nd), ep1(len, nd) !RomP
    8588  REAL epmlmmm1(len, nd, nd), eplamm1(len, nd) !RomP+jyg
     89  REAL qtc1(len, nd), sigt1(len, nd) !RomP
    8690  REAL wdtraina1(len, nd), wdtrainm1(len, nd) !RomP
    8791
    … …  
    141145      wdtraina1(idcum(i), k) = wdtraina(i, k) !RomP
    142146      wdtrainm1(idcum(i), k) = wdtrainm(i, k) !RomP
     147      qtc1(idcum(i), k) = qtc(i, k)
     148      sigt1(idcum(i), k) = sigt(i, k)
    143149
    144150    END DO

  • LMDZ5/branches/testing/libf/phylmd/cv3p1_closure.F90

    r1999 r2220  
    33
    44SUBROUTINE cv3p1_closure(nloc, ncum, nd, icb, inb, pbase, plcl, p, ph, tv, &
    5     tvp, buoy, supmax, ok_inhib, ale, alp, sig, w0, ptop2, cape, cin, m, &
     5    tvp, buoy, supmax, ok_inhib, ale, alp, omega,sig, w0, ptop2, cape, cin, m, &
    66    iflag, coef, plim1, plim2, asupmax, supmax0, asupmaxmin, cbmf, plfc, &
    77    wbeff)
    … …  
    3737  LOGICAL ok_inhib ! enable convection inhibition by dryness
    3838  REAL ale(nloc), alp(nloc)
     39  REAL omega(nloc,nd)
    3940
    4041  ! input/output:
    … …  
    5253
    5354  ! local variables:
    54   INTEGER il, i, j, k, icbmax, i0(nloc)
     55  INTEGER il, i, j, k, icbmax, i0(nloc), klfc
    5556  REAL deltap, fac, w, amu
    5657  REAL rhodp
    … …  
    523524  END DO
    524525
     526!CR:Compute k at plfc
     527  DO k=1,nl
     528     DO il=1,ncum
     529        if ((plfc(il).lt.ph(il,k)).and.(plfc(il).ge.ph(il,k+1))) then
     530           klfc=k
     531        endif
     532     ENDDO
     533  ENDDO
     534!RC
    525535
    526536  DO il = 1, ncum
    … …  
    528538    ! c       cbmf1(il) = alp2(il)/(0.5*wb*wb-Cin(il))
    529539    cbmf1(il) = alp2(il)/(2.*wbeff(il)*wbeff(il)-cin(il))
     540!CR: Add large-scale component to the mass-flux
     541!encore connu sous le nom "Experience du tube de dentifrice"
     542    if (coef_clos_ls.gt.0.) then
     543       cbmf1(il) = cbmf1(il) - coef_clos_ls*min(0.,1./RG*omega(il,klfc))
     544    endif
     545!RC
    530546    IF (cbmf1(il)==0 .AND. alp2(il)/=0.) THEN
    531547      WRITE (lunout, *) 'cv3p1_closure cbmf1=0 and alp NE 0 il alp2 alp cin ' &

  • LMDZ5/branches/testing/libf/phylmd/cva_driver.F90

    r2160 r2220  
    77                      u1, v1, tra1, &
    88                      p1, ph1, &
    9                       Ale1, Alp1, &
     9                      Ale1, Alp1, omega1, &
    1010                      sig1feed1, sig2feed1, wght1, &
    1111                      iflag1, ft1, fq1, fu1, fv1, ftra1, &
    … …  
    2424                      da1, phi1, mp1, phi21, d1a1, dam1, sigij1, wghti1, & ! RomP, RL
    2525                      clw1, elij1, evap1, ep1, epmlmMm1, eplaMm1, &        ! RomP, RL
    26                       wdtrainA1, wdtrainM1)                                ! RomP
     26                      wdtrainA1, wdtrainM1, qtc1, sigt1, tau_cld_cv, &
     27                      coefw_cld_cv)                                        ! RomP, AJ
    2728! **************************************************************
    2829! *
    … …  
    5556! iflag_ice_thermo Integer        Input        accounting for ice thermodynamics (0/1)
    5657! iflag_clos    Integer        Input        version of closure (0/1)
     58! tau_cld_cv    Real           Input        characteristic time of dissipation of mixing fluxes
     59! coefw_cld_cv  Real           Input        coefficient for updraft velocity in convection
    5760! ok_conserv_q  Logical        Input        when true corrections for water conservation are swtiched on
    5861! delt          Real           Input        time step
    … …  
    119122! phi1          Real           Output     used in tracer transport (cvltr)
    120123! mp1           Real           Output     used in tracer transport (cvltr)
    121                                          
     124! qtc1          Real           Output     specific humidity in convection
     125! sigt1         Real           Output     surface fraction in adiabatic updrafts                                         
    122126! phi21         Real           Output     used in tracer transport (cvltr)
    123127                                         
    … …  
    163167  INTEGER iflag_clos
    164168  LOGICAL ok_conserv_q
     169  REAL tau_cld_cv
     170  REAL coefw_cld_cv
    165171  REAL delt
    166172  REAL t1(len, nd)
    … …  
    178184  REAL Ale1(len)
    179185  REAL Alp1(len)
     186  REAL omega1(len,nd)
    180187  REAL sig1feed1 ! pressure at lower bound of feeding layer
    181188  REAL sig2feed1 ! pressure at upper bound of feeding layer
    … …  
    225232  REAL asupmaxmin1(len)
    226233  INTEGER lalim_conv(len)
     234  REAL qtc1(len, nd)         ! cld
     235  REAL sigt1(len, nd)        ! cld
     236
    227237! RomP >>>
    228238  REAL wdtrainA1(len, nd), wdtrainM1(len, nd)
    … …  
    455465  REAL supmax(nloc, klev)
    456466  REAL Ale(nloc), Alp(nloc), coef_clos(nloc)
     467  REAL omega(nloc,klev)
    457468  REAL sigd(nloc)
    458469! real mp(nloc,klev), qp(nloc,klev), up(nloc,klev), vp(nloc,klev)
    … …  
    487498  REAL hnk(nloc), unk(nloc), vnk(nloc)
    488499
     500  REAL qtc(nloc, klev)         ! cld
     501  REAL sigt(nloc, klev)        ! cld
     502 
    489503! RomP >>>
    490504  REAL wdtrainA(nloc, klev), wdtrainM(nloc, klev)
    … …  
    593607
    594608! RomP >>>
     609  sigt1(:, :) = 0.
     610  qtc1(:, :) = 0.
    595611  wdtrainA1(:, :) = 0.
    596612  wdtrainM1(:, :) = 0.
    … …  
    776792                         h1_wake, lv1_wake, lf1_wake, cpn1_wake, tv1_wake, &
    777793                         sig1, w01, ptop21, &
    778                          Ale1, Alp1, &
     794                         Ale1, Alp1, omega1, &
    779795                         iflag, nk, icb, icbs, &
    780796                         plcl, tnk, qnk, gznk, hnk, unk, vnk, &
    … …  
    786802                         h_wake, lv_wake, lf_wake, cpn_wake, tv_wake, &
    787803                         sig, w0, ptop2, &
    788                          Ale, Alp)
     804                         Ale, Alp, omega)
    789805
    790806! print*,'tv ',tv
    … …  
    877893        CALL cv3p1_closure(nloc, ncum, nd, icb, inb, &         ! na->nd
    878894                           pbase, plcl, p, ph, tv, tvp, buoy, &
    879                            supmax, ok_inhib, Ale, Alp, &
     895                           supmax, ok_inhib, Ale, Alp, omega, &
    880896                           sig, w0, ptop2, cape, cin, m, iflag, coef_clos, &
    881897                           Plim1, plim2, asupmax, supmax0, &
    … …  
    978994                     cbmf, upwd, dnwd, dnwd0, ma, mip, &
    979995                     tls, tps, qcondc, wd, &
    980                      ftd, fqd)
     996                     ftd, fqd, qnk, qtc, sigt, tau_cld_cv, coefw_cld_cv)
    981997    END IF
    982998
    … …  
    10321048                           clw, elij, evap, ep, epmlmMm, eplaMm, &       ! RomP
    10331049                           wdtrainA, wdtrainM, &                         ! RomP
     1050                           qtc, sigt, &
    10341051                           iflag1, kbas1, ktop1, &
    10351052                           precip1, cbmf1, plcl1, plfc1, wbeff1, sig1, w01, ptop21, &
    … …  
    10431060                           da1, phi1, mp1, phi21, d1a1, dam1, sigij1,  & ! RomP
    10441061                           clw1, elij1, evap1, ep1, epmlmMm1, eplaMm1, & ! RomP
    1045                            wdtrainA1, wdtrainM1)                         ! RomP
     1062                           wdtrainA1, wdtrainM1,                       & ! RomP
     1063                           qtc1, sigt1)
    10461064    END IF
    10471065

  • LMDZ5/branches/testing/libf/phylmd/fisrtilp.F90

    r2160 r2220  
    88     frac_impa, frac_nucl, beta,                        &
    99     prfl, psfl, rhcl, zqta, fraca,                     &
    10      ztv, zpspsk, ztla, zthl, iflag_cldcon,             &
     10     ztv, zpspsk, ztla, zthl, iflag_cldth,             &
    1111     iflag_ice_thermo)
    1212
    … …  
    8282  INTEGER ninter ! sous-intervals pour la precipitation
    8383  INTEGER ncoreczq 
    84   INTEGER iflag_cldcon
     84  INTEGER iflag_cldth
    8585  INTEGER iflag_ice_thermo
    8686  PARAMETER (ninter=5)
    … …  
    545545           enddo
    546546
    547            if (iflag_cldcon>=5) then
     547           if (iflag_cldth>=5) then
    548548
    549549              call cloudth(klon,klev,k,ztv, &
    … …  
    559559           endif
    560560
    561            if (iflag_cldcon <= 4) then
     561           if (iflag_cldth <= 4) then
    562562              lognormale = .true.
    563            elseif (iflag_cldcon >= 6) then
     563           elseif (iflag_cldth >= 6) then
    564564              ! lognormale en l'absence des thermiques
    565565              lognormale = fraca(:,k) < 1e-10
    566566           else
    567               ! Dans le cas iflag_cldcon=5, on prend systématiquement la
     567              ! Dans le cas iflag_cldth=5, on prend systématiquement la
    568568              ! bi-gaussienne
    569569              lognormale = .false.

  • LMDZ5/branches/testing/libf/phylmd/hines_gwd.F90

    r1999 r2220  
    625625    mmin_alpha, kstar, slope, f1, f2, f3, naz, levbot, levtop, il1, il2, &
    626626    nlons, nlevs, nazmth, sigsqmcw, sigmatm, lorms, sigalpmc, f2mod)
    627 
    628627  ! Smooth cutoff wavenumbers and total rms velocity in the vertical
    629628  ! direction NSMAX times, using FLUX_U as temporary work array.
    … …  
    715714    i_alpha, mmin_alpha, kstar, slope, f1, f2, f3, naz, levbot, levtop, il1, &
    716715    il2, nlons, nlevs, nazmth, sigsqmcw, sigmatm, lorms, sigalpmc, f2mod)
    717 
     716  IMPLICIT NONE
    718717  ! This routine calculates the cutoff vertical wavenumber and velocity
    719718  ! variances on a longitude by altitude grid for the Hines' Doppler
    … …  
    766765
    767766  INTEGER naz, levbot, levtop, il1, il2, nlons, nlevs, nazmth
    768   REAL slope, kstar(nlons), f1, f2, f3
     767  REAL slope, kstar(nlons), f1, f2, f3, f2mfac
    769768  REAL m_alpha(nlons, nlevs, nazmth)
    770769  REAL sigma_alpha(nlons, nlevs, nazmth)
    … …  
    938937SUBROUTINE hines_wind(v_alpha, vel_u, vel_v, naz, il1, il2, lev1, lev2, &
    939938    nlons, nlevs, nazmth)
    940 
     939    IMPLICIT NONE
    941940  ! This routine calculates the azimuthal horizontal background wind
    942941  ! components
    … …  
    10341033    m_alpha, ak_alpha, k_alpha, slope, naz, il1, il2, lev1, lev2, nlons, &
    10351034    nlevs, nazmth, lorms)
    1036 
     1035    IMPLICIT NONE
    10371036  ! Calculate zonal and meridional components of the vertical flux
    10381037  ! of horizontal momentum and corresponding wave drag (force per unit mass)
    … …  
    10891088  ! Internal variables.
    10901089
    1091   INTEGER i, l, lev1p, lev2m
     1090  INTEGER i, l, lev1p, lev2m, lev2p
    10921091  REAL cos45, prod2, prod4, prod6, prod8, dendz, dendz2
    10931092  DATA cos45/0.7071068/
    … …  
    12341233    bvfreq, density, densb, sigma_t, visc_mol, kstar, slope, f2, f3, f5, f6, &
    12351234    naz, il1, il2, lev1, lev2, nlons, nlevs, nazmth)
    1236 
     1235  IMPLICIT NONE
    12371236  ! This routine calculates the gravity wave induced heating and
    12381237  ! diffusion coefficient on a longitude by altitude grid for
    … …  
    13551354SUBROUTINE hines_sigma(sigma_t, sigma_alpha, sigsqh_alpha, naz, lev, il1, &
    13561355    il2, nlons, nlevs, nazmth)
    1357 
     1356  IMPLICIT NONE
    13581357  ! This routine calculates the total rms and azimuthal rms horizontal
    13591358  ! velocities at a given level on a longitude by altitude grid for
    … …  
    14501449SUBROUTINE hines_intgrl(i_alpha, v_alpha, m_alpha, bvfb, slope, naz, lev, &
    14511450    il1, il2, nlons, nlevs, nazmth, lorms)
    1452 
     1451  IMPLICIT NONE
    14531452  ! This routine calculates the vertical wavenumber integral
    14541453  ! for a single vertical level at each azimuth on a longitude grid
    … …  
    16431642    alt_cutoff, smco, nsmax, iheatcal, k_alpha, ierror, nmessg, nlons, &
    16441643    nazmth, coslat)
    1645 
     1644  IMPLICIT NONE
    16461645  ! This routine specifies various parameters needed for the
    16471646  ! the Hines' Doppler spread gravity wave drag parameterization scheme.
    … …  
    17721771    sigma_alpha, v_alpha, m_alpha, iu_print, iv_print, nmessg, ilprt1, &
    17731772    ilprt2, levprt1, levprt2, naz, nlons, nlevs, nazmth)
    1774 
     1773  IMPLICIT NONE
    17751774  ! Print out altitude profiles of various quantities from
    17761775  ! Hines' Doppler spread gravity wave drag parameterization scheme.
    … …  
    18641863SUBROUTINE hines_exp(data, data_zmax, alt, alt_exp, iorder, il1, il2, lev1, &
    18651864    lev2, nlons, nlevs)
    1866 
     1865  IMPLICIT NONE
    18671866  ! This routine exponentially damps a longitude by altitude array
    18681867  ! of data above a specified altitude.
    … …  
    19411940SUBROUTINE vert_smooth(data, work, coeff, nsmooth, il1, il2, lev1, lev2, &
    19421941    nlons, nlevs)
    1943 
     1942  IMPLICIT NONE
    19441943  ! Smooth a longitude by altitude array in the vertical over a
    19451944  ! specified number of levels using a three point smoother.

  • LMDZ5/branches/testing/libf/phylmd/ini_wake.F90

    r1999 r2220  
    44SUBROUTINE ini_wake(wape, fip, it_wape_prescr, wape_prescr, fip_prescr, &
    55    alp_bl_prescr, ale_bl_prescr)
     6  IMPLICIT NONE
    67  ! **************************************************************
    78  ! *
    … …  
    3940  include 'iniprint.h'
    4041  ! declarations
     42  REAL wape, fip, wape_prescr, fip_prescr
     43  INTEGER it_wape_prescr
    4144  REAL ale_bl_prescr
    4245  REAL alp_bl_prescr
    4346  REAL it
     47  REAL w,f,alebl,alpbl
    4448
    4549  ! FH A mettre si besoin dans physiq.def

  • LMDZ5/branches/testing/libf/phylmd/limit_read_mod.F90

    r1910 r2220  
    9393!****************************************************************************************
    9494
    95     IF (type_ocean == 'couple') THEN
     95IF (type_ocean == 'couple'.OR. &
     96         (type_ocean == 'slab' .AND. version_ocean == 'sicINT')) THEN
    9697       ! limit.nc has not yet been read. Do it now!
    9798       CALL limit_read_tot(itime, dtime, jour, is_modified)

  • LMDZ5/branches/testing/libf/phylmd/limit_slab.F90

    r2073 r2220  
    11! $Header$
    22
    3 SUBROUTINE limit_slab(itime, dtime, jour, lmt_bils, diff_sst)
     3SUBROUTINE limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv)
    44
    55  USE dimphy
    … …  
    2020  INTEGER, INTENT(IN) :: jour    ! jour a lire dans l'annee
    2121  REAL   , INTENT(IN) :: dtime   ! pas de temps de la physique (en s)
    22   REAL, DIMENSION(klon), INTENT(OUT) :: lmt_bils, diff_sst
     22  REAL, DIMENSION(klon), INTENT(OUT) :: lmt_bils, diff_sst, diff_siv
    2323
    2424! Locals variables with attribute SAVE
    2525!****************************************************************************************
    2626  REAL, DIMENSION(:), ALLOCATABLE, SAVE :: bils_save, diff_sst_save
    27 !$OMP THREADPRIVATE(bils_save, diff_sst_save)
     27  REAL, DIMENSION(:), ALLOCATABLE, SAVE :: diff_siv_save
     28!$OMP THREADPRIVATE(bils_save, diff_sst_save, diff_siv_save)
    2829
    2930! Locals variables
    … …  
    3334  INTEGER, DIMENSION(2)    :: start, epais
    3435  REAL, DIMENSION(klon_glo):: bils_glo, sst_l_glo, sst_lp1_glo, diff_sst_glo
     36  REAL, DIMENSION(klon_glo):: siv_l_glo, siv_lp1_glo, diff_siv_glo
    3537  CHARACTER (len = 20)     :: modname = 'limit_slab'
    36   LOGICAL                  :: read_bils,read_sst
     38  LOGICAL                  :: read_bils,read_sst,read_siv
    3739
    3840! End declaration
    … …  
    4951        read_bils=.TRUE.
    5052        read_sst=.TRUE.
     53        read_siv=.TRUE.
    5154       
    5255        ierr = NF90_OPEN ('limit_slab.nc', NF90_NOWRITE, nid)
    … …  
    5457            read_bils=.FALSE.
    5558            read_sst=.FALSE.
     59            read_siv=.FALSE.
    5660        ELSE ! read file
    5761       
    … …  
    6367
    6468!****************************************************************************************
    65 ! 2) Read bils and SST tendency
     69! 2) Read bils and SST/ ice volume tendency
    6670!
    6771!****************************************************************************************
    … …  
    8993        END IF
    9094
     95! Read siv_glo for this day
     96        ierr = NF90_INQ_VARID(nid, 'SICV', nvarid)
     97        IF (ierr /= NF90_NOERR)  THEN
     98            read_siv=.FALSE.
     99        ELSE
     100            start(2) = jour
     101            ierr = NF90_GET_VAR(nid,nvarid,siv_l_glo,start,epais)
     102            IF (ierr /= NF90_NOERR) read_siv=.FALSE.
     103! Read siv_glo for one day ahead
     104            start(2) = jour + 1
     105            IF (start(2) > 360) start(2)=1
     106            ierr = NF90_GET_VAR(nid,nvarid,siv_lp1_glo,start,epais)
     107            IF (ierr /= NF90_NOERR) read_siv=.FALSE.
     108        END IF
     109
    91110!****************************************************************************************
    92111! 5) Close file and distribute variables to all processus
    … …  
    98117        IF (read_sst) THEN
    99118! Calculate difference in temperature between this day and one ahead
    100 !            DO i=1, klon_glo-1
    101 !               diff_sst_glo(i) = sst_lp1_glo(i) - sst_l_glo(i)
    102 !            END DO
    103 !            diff_sst_glo(klon_glo) = sst_lp1_glo(klon_glo) - sst_l_glo(1)
    104119            DO i=1, klon_glo
    105120               diff_sst_glo(i) = sst_lp1_glo(i) - sst_l_glo(i)
    106121            END DO
    107122        END IF !read_sst
     123        IF (read_siv) THEN
     124! Calculate difference in temperature between this day and one ahead
     125            DO i=1, klon_glo
     126               diff_siv_glo(i) = siv_lp1_glo(i) - siv_l_glo(i)
     127            END DO
     128        END IF !read_siv
    108129     ENDIF ! is_mpi_root
    109130
    … …  
    111132       
    112133     IF (.NOT. ALLOCATED(bils_save)) THEN
    113         ALLOCATE(bils_save(klon), diff_sst_save(klon), stat=ierr)
     134        ALLOCATE(bils_save(klon), diff_sst_save(klon), diff_siv_save(klon), stat=ierr)
    114135        IF (ierr /= 0) CALL abort_gcm('limit_slab', 'pb in allocation',1)
    115136     END IF
    116137
    117 ! Giveddefault values if needed
     138! Give default values if needed
    118139     IF (read_bils) THEN
    119140         CALL Scatter(bils_glo, bils_save)
    … …  
    126147         diff_sst_save(:)=0.
    127148     END IF
     149     IF (read_siv) THEN
     150         CALL Scatter(diff_siv_glo, diff_siv_save)
     151     ELSE
     152         diff_siv_save(:)=0.
     153     END IF
    128154     
    129155  ENDIF ! time to read
    … …  
    131157  lmt_bils(:) = bils_save(:)
    132158  diff_sst(:) = diff_sst_save(:)
     159  diff_siv(:) = diff_siv_save(:)
     160
    133161 
    134162END SUBROUTINE limit_slab

  • LMDZ5/branches/testing/libf/phylmd/lmdz1d.F90

    r2187 r2220  
    2222      USE phyaqua_mod
    2323      USE mod_1D_cases_read
     24      USE mod_1D_amma_read
    2425
    2526      implicit none
    … …  
    118119        logical :: forcing_astex   = .false.
    119120        logical :: forcing_fire    = .false.
     121        logical :: forcing_case    = .false.
    120122        integer :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file
    121123!                                                            (cf read_tsurf1d.F)
    … …  
    171173      real :: dt_cooling(llm),d_th_adv(llm),d_t_nudge(llm)
    172174      real :: d_u_nudge(llm),d_v_nudge(llm)
     175      real :: du_adv(llm),dv_adv(llm)
    173176      real :: alpha
    174177      real :: ttt
    … …  
    285288!             Different stages: soil model alone, atm. model alone
    286289!             then both models coupled
     290!forcing_type = 10 ==> forcing_case = .true.
     291!             initial profiles and large scale forcings in cas.nc
     292!             LS convergence, omega and SST imposed from CINDY-DYNAMO files
    287293!forcing_type = 40 ==> forcing_GCSSold = .true.
    288294!             initial profile from GCSS file
    … …  
    321327      elseif (forcing_type .eq.7) THEN
    322328       forcing_dice = .true.
     329      elseif (forcing_type .eq.10) THEN
     330       forcing_case = .true.
    323331      elseif (forcing_type .eq.40) THEN
    324332       forcing_GCSSold = .true.
    … …  
    428436     & (year_ini_dice,mth_ini_dice,day_ini_dice,heure_ini_dice             &
    429437     & ,day_ju_ini_dice)
     438      ELSEIF (forcing_type .eq.10) THEN
     439! Convert the initial date to Julian day
     440      print*,'time cindy',year_ini_cas,mth_ini_cas,day_ini_cas
     441      call ymds2ju                                                         &
     442     & (year_ini_cas,mth_ini_cas,day_ini_cas,heure_ini_cas &
     443     & ,day_ju_ini_cas)
     444      print*,'time cindy 2',day_ju_ini_cas
    430445      ELSEIF (forcing_type .eq.59) THEN
    431446! Convert the initial date of Sandu case to Julian day
    … …  
    943958        u(1:mxcalc)=u(1:mxcalc) + timestep*(                                &
    944959     &              du_phys(1:mxcalc)                                       &
    945      &             +du_age(1:mxcalc)                                        &
     960     &             +du_age(1:mxcalc)+du_adv(1:mxcalc)                       &
    946961     &             +d_u_nudge(1:mxcalc) )           
    947962        v(1:mxcalc)=v(1:mxcalc) + timestep*(                                 &
    948963     &              dv_phys(1:mxcalc)                                       &
    949      &             +dv_age(1:mxcalc)                                        &
     964     &             +dv_age(1:mxcalc)+dv_adv(1:mxcalc)                       &
    950965     &             +d_v_nudge(1:mxcalc) )
    951966        q(1:mxcalc,:)=q(1:mxcalc,:)+timestep*(                              &

  • LMDZ5/branches/testing/libf/phylmd/mod_1D_cases_read.F90

    r2160 r2220  
    22
    33!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    4 !Declarations specifiques au cas AMMA
    5         character*80 :: fich_amma
    6 ! Option du cas AMMA ou on impose la discretisation verticale (Ap,Bp)
    7         integer nlev_amma, nt_amma
    8 
    9 
    10         integer year_ini_amma, day_ini_amma, mth_ini_amma
    11         real heure_ini_amma
    12         real day_ju_ini_amma   ! Julian day of amma first day
    13         parameter (year_ini_amma=2006)
    14         parameter (mth_ini_amma=7)
    15         parameter (day_ini_amma=10)  ! 10 = 10Juil2006
    16         parameter (heure_ini_amma=0.) !0h en secondes
    17         real dt_amma
    18         parameter (dt_amma=1800.)
    19 
    20 !profils initiaux:
    21         real, allocatable::  plev_amma(:)
    22 
    23         real, allocatable::  z_amma(:)
    24         real, allocatable::  th_amma(:),q_amma(:)
    25         real, allocatable::  u_amma(:)
    26         real, allocatable::  v_amma(:)
    27 
    28         real, allocatable::  th_ammai(:),q_ammai(:)
    29         real, allocatable::  u_ammai(:)
    30         real, allocatable::  v_ammai(:)
    31         real, allocatable::  vitw_ammai(:)
    32         real, allocatable::  ht_ammai(:)
    33         real, allocatable::  hq_ammai(:)
    34         real, allocatable::  vt_ammai(:)
    35         real, allocatable::  vq_ammai(:)
    36 
    37 !forcings
    38         real, allocatable::  ht_amma(:,:)
    39         real, allocatable::  hq_amma(:,:)
    40         real, allocatable::  vitw_amma(:,:)
    41         real, allocatable::  lat_amma(:),sens_amma(:)
     4!Declarations specifiques au cas standard
     5        character*80 :: fich_cas
     6! Discrétisation
     7        integer nlev_cas, nt_cas
     8
     9
     10        integer year_ini_cas, day_ini_cas, mth_ini_cas
     11        real heure_ini_cas
     12        real day_ju_ini_cas   ! Julian day of case first day
     13        parameter (year_ini_cas=2011)
     14        parameter (mth_ini_cas=10)
     15        parameter (day_ini_cas=1)  ! 10 = 10Juil2006
     16        parameter (heure_ini_cas=0.) !0h en secondes
     17        real pdt_cas
     18        parameter (pdt_cas=3.*3600)
     19
     20!CR ATTENTION TEST AMMA
     21!        parameter (year_ini_cas=2006)
     22!        parameter (mth_ini_cas=7)
     23!        parameter (day_ini_cas=10)  ! 10 = 10Juil2006
     24!        parameter (heure_ini_cas=0.) !0h en secondes
     25!        parameter (pdt_cas=1800.)
     26
     27!profils environnementaux
     28        real, allocatable::  plev_cas(:,:)
     29
     30        real, allocatable::  z_cas(:,:)
     31        real, allocatable::  t_cas(:,:),q_cas(:,:),rh_cas(:,:)
     32        real, allocatable::  th_cas(:,:),rv_cas(:,:)
     33        real, allocatable::  u_cas(:,:)
     34        real, allocatable::  v_cas(:,:)
     35
     36!forcing
     37        real, allocatable::  ht_cas(:,:),vt_cas(:,:),dt_cas(:,:),dtrad_cas(:,:)
     38        real, allocatable::  hth_cas(:,:),vth_cas(:,:),dth_cas(:,:)
     39        real, allocatable::  hq_cas(:,:),vq_cas(:,:),dq_cas(:,:)
     40        real, allocatable::  hr_cas(:,:),vr_cas(:,:),dr_cas(:,:)
     41        real, allocatable::  hu_cas(:,:),vu_cas(:,:),du_cas(:,:)
     42        real, allocatable::  hv_cas(:,:),vv_cas(:,:),dv_cas(:,:)
     43        real, allocatable::  vitw_cas(:,:)
     44        real, allocatable::  ug_cas(:,:),vg_cas(:,:)
     45        real, allocatable::  lat_cas(:),sens_cas(:),ts_cas(:)
    4246
    4347!champs interpoles
    44         real, allocatable::  vitw_profamma(:)
    45         real, allocatable::  ht_profamma(:)
    46         real, allocatable::  hq_profamma(:)
    47         real lat_profamma,sens_profamma
    48         real, allocatable::  vt_profamma(:)
    49         real, allocatable::  vq_profamma(:)
    50         real, allocatable::  th_profamma(:)
    51         real, allocatable::  q_profamma(:)
    52         real, allocatable::  u_profamma(:)
    53         real, allocatable::  v_profamma(:)
     48        real, allocatable::  plev_prof_cas(:)
     49        real, allocatable::  t_prof_cas(:)
     50        real, allocatable::  q_prof_cas(:)
     51        real, allocatable::  u_prof_cas(:)
     52        real, allocatable::  v_prof_cas(:)       
     53
     54        real, allocatable::  vitw_prof_cas(:)
     55        real, allocatable::  ug_prof_cas(:)
     56        real, allocatable::  vg_prof_cas(:)
     57        real, allocatable::  ht_prof_cas(:)
     58        real, allocatable::  hq_prof_cas(:)
     59        real, allocatable::  vt_prof_cas(:)
     60        real, allocatable::  vq_prof_cas(:)
     61        real, allocatable::  dt_prof_cas(:)
     62        real, allocatable::  dtrad_prof_cas(:)
     63        real, allocatable::  dq_prof_cas(:)
     64        real, allocatable::  hu_prof_cas(:)
     65        real, allocatable::  hv_prof_cas(:)
     66        real, allocatable::  vu_prof_cas(:)
     67        real, allocatable::  vv_prof_cas(:)
     68        real, allocatable::  du_prof_cas(:)
     69        real, allocatable::  dv_prof_cas(:)
     70
     71
     72        real lat_prof_cas,sens_prof_cas,ts_prof_cas
     73     
    5474
    5575
    5676CONTAINS
    5777
    58 SUBROUTINE read_1D_cases
     78SUBROUTINE read_1D_cas
    5979      implicit none
    6080
    … …  
    6282
    6383      INTEGER nid,rid,ierr
    64 
    65       fich_amma='amma.nc'
    66       print*,'fich_amma ',fich_amma
    67       ierr = NF_OPEN(fich_amma,NF_NOWRITE,nid)
    68       print*,'fich_amma,NF_NOWRITE,nid ',fich_amma,NF_NOWRITE,nid
     84      INTEGER ii,jj
     85
     86      fich_cas='setup/cas.nc'
     87      print*,'fich_cas ',fich_cas
     88      ierr = NF_OPEN(fich_cas,NF_NOWRITE,nid)
     89      print*,'fich_cas,NF_NOWRITE,nid ',fich_cas,NF_NOWRITE,nid
    6990      if (ierr.NE.NF_NOERR) then
    7091         write(*,*) 'ERROR: GROS Pb opening forcings nc file '
    … …  
    7394      endif
    7495!.......................................................................
     96      ierr=NF_INQ_DIMID(nid,'lat',rid)
     97      IF (ierr.NE.NF_NOERR) THEN
     98         print*, 'Oh probleme lecture dimension lat'
     99      ENDIF
     100      ierr=NF_INQ_DIMLEN(nid,rid,ii)
     101      print*,'OK nid,rid,lat',nid,rid,ii
     102!.......................................................................
     103      ierr=NF_INQ_DIMID(nid,'lon',rid)
     104      IF (ierr.NE.NF_NOERR) THEN
     105         print*, 'Oh probleme lecture dimension lon'
     106      ENDIF
     107      ierr=NF_INQ_DIMLEN(nid,rid,jj)
     108      print*,'OK nid,rid,lat',nid,rid,jj
     109!.......................................................................
    75110      ierr=NF_INQ_DIMID(nid,'lev',rid)
    76111      IF (ierr.NE.NF_NOERR) THEN
    77112         print*, 'Oh probleme lecture dimension zz'
    78113      ENDIF
    79       ierr=NF_INQ_DIMLEN(nid,rid,nlev_amma)
    80       print*,'OK nid,rid,nlev_amma',nid,rid,nlev_amma
     114      ierr=NF_INQ_DIMLEN(nid,rid,nlev_cas)
     115      print*,'OK nid,rid,nlev_cas',nid,rid,nlev_cas
    81116!.......................................................................
    82117      ierr=NF_INQ_DIMID(nid,'time',rid)
    83118      print*,'nid,rid',nid,rid
    84       nt_amma=0
     119      nt_cas=0
    85120      IF (ierr.NE.NF_NOERR) THEN
    86121        stop 'probleme lecture dimension sens'
    87122      ENDIF
    88       ierr=NF_INQ_DIMLEN(nid,rid,nt_amma)
    89       print*,'nid,rid,nlev_amma',nid,rid,nt_amma
     123      ierr=NF_INQ_DIMLEN(nid,rid,nt_cas)
     124      print*,'nid,rid,nlev_cas',nid,rid,nt_cas
    90125
    91126!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    92 !profils initiaux:
    93         allocate(plev_amma(nlev_amma))
    94        
    95         allocate(z_amma(nlev_amma))
    96         allocate(th_amma(nlev_amma),q_amma(nlev_amma))
    97         allocate(u_amma(nlev_amma))
    98         allocate(v_amma(nlev_amma))
    99 
    100 !forcings
    101         allocate(ht_amma(nlev_amma,nt_amma))
    102         allocate(hq_amma(nlev_amma,nt_amma))
    103         allocate(vitw_amma(nlev_amma,nt_amma))
    104         allocate(lat_amma(nt_amma),sens_amma(nt_amma))
    105 
    106 !profils initiaux:
    107         allocate(th_ammai(nlev_amma),q_ammai(nlev_amma))
    108         allocate(u_ammai(nlev_amma))
    109         allocate(v_ammai(nlev_amma))
    110         allocate(vitw_ammai(nlev_amma) )
    111         allocate(ht_ammai(nlev_amma))
    112         allocate(hq_ammai(nlev_amma))
    113         allocate(vt_ammai(nlev_amma))
    114         allocate(vq_ammai(nlev_amma))
     127!profils moyens:
     128        allocate(plev_cas(nlev_cas,nt_cas))       
     129        allocate(z_cas(nlev_cas,nt_cas))
     130        allocate(t_cas(nlev_cas,nt_cas),q_cas(nlev_cas,nt_cas),rh_cas(nlev_cas,nt_cas))
     131        allocate(th_cas(nlev_cas,nt_cas),rv_cas(nlev_cas,nt_cas))
     132        allocate(u_cas(nlev_cas,nt_cas))
     133        allocate(v_cas(nlev_cas,nt_cas))
     134
     135!forcing
     136        allocate(ht_cas(nlev_cas,nt_cas),vt_cas(nlev_cas,nt_cas),dt_cas(nlev_cas,nt_cas),dtrad_cas(nlev_cas,nt_cas))
     137        allocate(hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas),dq_cas(nlev_cas,nt_cas))
     138        allocate(hth_cas(nlev_cas,nt_cas),vth_cas(nlev_cas,nt_cas),dth_cas(nlev_cas,nt_cas))
     139        allocate(hr_cas(nlev_cas,nt_cas),vr_cas(nlev_cas,nt_cas),dr_cas(nlev_cas,nt_cas))
     140        allocate(hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas),du_cas(nlev_cas,nt_cas))
     141        allocate(hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas),dv_cas(nlev_cas,nt_cas))
     142        allocate(vitw_cas(nlev_cas,nt_cas))
     143        allocate(ug_cas(nlev_cas,nt_cas))
     144        allocate(vg_cas(nlev_cas,nt_cas))
     145        allocate(lat_cas(nt_cas),sens_cas(nt_cas),ts_cas(nt_cas))
     146
    115147
    116148!champs interpoles
    117         allocate(vitw_profamma(nlev_amma))
    118         allocate(ht_profamma(nlev_amma))
    119         allocate(hq_profamma(nlev_amma))
    120         allocate(vt_profamma(nlev_amma))
    121         allocate(vq_profamma(nlev_amma))
    122         allocate(th_profamma(nlev_amma))
    123         allocate(q_profamma(nlev_amma))
    124         allocate(u_profamma(nlev_amma))
    125         allocate(v_profamma(nlev_amma))
     149        allocate(plev_prof_cas(nlev_cas))
     150        allocate(t_prof_cas(nlev_cas))
     151        allocate(q_prof_cas(nlev_cas))
     152        allocate(u_prof_cas(nlev_cas))
     153        allocate(v_prof_cas(nlev_cas))
     154
     155        allocate(vitw_prof_cas(nlev_cas))
     156        allocate(ug_prof_cas(nlev_cas))
     157        allocate(vg_prof_cas(nlev_cas))
     158        allocate(ht_prof_cas(nlev_cas))
     159        allocate(hq_prof_cas(nlev_cas))
     160        allocate(hu_prof_cas(nlev_cas))
     161        allocate(hv_prof_cas(nlev_cas))
     162        allocate(vt_prof_cas(nlev_cas))
     163        allocate(vq_prof_cas(nlev_cas))
     164        allocate(vu_prof_cas(nlev_cas))
     165        allocate(vv_prof_cas(nlev_cas))
     166        allocate(dt_prof_cas(nlev_cas))
     167        allocate(dtrad_prof_cas(nlev_cas))
     168        allocate(dq_prof_cas(nlev_cas))
     169        allocate(du_prof_cas(nlev_cas))
     170        allocate(dv_prof_cas(nlev_cas))
    126171
    127172        print*,'Allocations OK'
    128         call read_amma(nid,nlev_amma,nt_amma                                  &
    129      &     ,z_amma,plev_amma,th_amma,q_amma,u_amma,v_amma,vitw_amma         &
    130      &     ,ht_amma,hq_amma,sens_amma,lat_amma)
    131 
    132 END SUBROUTINE read_1D_cases
     173        call read_cas(nid,nlev_cas,nt_cas                                  &
     174     &     ,z_cas,plev_cas,t_cas,q_cas,rh_cas,th_cas,rv_cas,u_cas,v_cas     &
     175     &     ,ug_cas,vg_cas,vitw_cas,du_cas,hu_cas,vu_cas,dv_cas,hv_cas,vv_cas &
     176     &     ,dt_cas,dtrad_cas,ht_cas,vt_cas,dq_cas,hq_cas,vq_cas             &
     177     &     ,dth_cas,hth_cas,vth_cas,dr_cas,hr_cas,vr_cas,sens_cas,lat_cas,ts_cas)
     178
     179
     180END SUBROUTINE read_1D_cas
    133181
    134182
    … …  
    136184!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    137185SUBROUTINE deallocate_1D_cases
    138 !profils initiaux:
    139         deallocate(plev_amma)
     186!profils environnementaux:
     187        deallocate(plev_cas)
    140188       
    141         deallocate(z_amma)
    142         deallocate(th_amma,q_amma)
    143         deallocate(u_amma)
    144         deallocate(v_amma)
    145 
    146         deallocate(th_ammai,q_ammai)
    147         deallocate(u_ammai)
    148         deallocate(v_ammai)
    149         deallocate(vitw_ammai )
    150         deallocate(ht_ammai)
    151         deallocate(hq_ammai)
    152         deallocate(vt_ammai)
    153         deallocate(vq_ammai)
     189        deallocate(z_cas)
     190        deallocate(t_cas,q_cas,rh_cas)
     191        deallocate(th_cas,rv_cas)
     192        deallocate(u_cas)
     193        deallocate(v_cas)
    154194       
    155 !forcings
    156         deallocate(ht_amma)
    157         deallocate(hq_amma)
    158         deallocate(vitw_amma)
    159         deallocate(lat_amma,sens_amma)
     195!forcing
     196        deallocate(ht_cas,vt_cas,dt_cas,dtrad_cas)
     197        deallocate(hq_cas,vq_cas,dq_cas)
     198        deallocate(hth_cas,vth_cas,dth_cas)
     199        deallocate(hr_cas,vr_cas,dr_cas)
     200        deallocate(hu_cas,vu_cas,du_cas)
     201        deallocate(hv_cas,vv_cas,dv_cas)
     202        deallocate(vitw_cas)
     203        deallocate(ug_cas)
     204        deallocate(vg_cas)
     205        deallocate(lat_cas,sens_cas,ts_cas)
    160206
    161207!champs interpoles
    162         deallocate(vitw_profamma)
    163         deallocate(ht_profamma)
    164         deallocate(hq_profamma)
    165         deallocate(vt_profamma)
    166         deallocate(vq_profamma)
    167         deallocate(th_profamma)
    168         deallocate(q_profamma)
    169         deallocate(u_profamma)
    170         deallocate(v_profamma)
     208        deallocate(plev_prof_cas)
     209        deallocate(t_prof_cas)
     210        deallocate(q_prof_cas)
     211        deallocate(u_prof_cas)
     212        deallocate(v_prof_cas)
     213
     214        deallocate(vitw_prof_cas)
     215        deallocate(ug_prof_cas)
     216        deallocate(vg_prof_cas)
     217        deallocate(ht_prof_cas)
     218        deallocate(hq_prof_cas)
     219        deallocate(hu_prof_cas)
     220        deallocate(hv_prof_cas)
     221        deallocate(vt_prof_cas)
     222        deallocate(vq_prof_cas)
     223        deallocate(vu_prof_cas)
     224        deallocate(vv_prof_cas)
     225        deallocate(dt_prof_cas)
     226        deallocate(dtrad_prof_cas)
     227        deallocate(dq_prof_cas)
     228        deallocate(du_prof_cas)
     229        deallocate(dv_prof_cas)
     230        deallocate(t_prof_cas)
     231        deallocate(q_prof_cas)
     232        deallocate(u_prof_cas)
     233        deallocate(v_prof_cas)
     234
    171235END SUBROUTINE deallocate_1D_cases
    172236
    … …  
    174238END MODULE mod_1D_cases_read
    175239!=====================================================================
    176       subroutine read_amma(nid,nlevel,ntime                          &
    177      &     ,zz,pp,temp,qv,u,v,dw                   &
    178      &     ,dt,dq,sens,flat)
    179 
    180 !program reading forcings of the AMMA case study
     240      subroutine read_cas(nid,nlevel,ntime                          &
     241     &     ,zz,pp,temp,qv,rh,theta,rv,u,v,ug,vg,w,                   &
     242     &     du,hu,vu,dv,hv,vv,dt,dtrad,ht,vt,dq,hq,vq,                     &
     243     &     dth,hth,vth,dr,hr,vr,sens,flat,ts)
     244
     245!program reading forcing of the case study
    181246      implicit none
    182247#include "netcdf.inc"
    … …  
    184249      integer ntime,nlevel
    185250
    186       real zz(nlevel)
    187       real temp(nlevel),pp(nlevel)
    188       real qv(nlevel),u(nlevel)
    189       real v(nlevel)
    190       real dw(nlevel,ntime)
    191       real dt(nlevel,ntime)
    192       real dq(nlevel,ntime)
    193       real flat(ntime),sens(ntime)
     251      real zz(nlevel,ntime)
     252      real pp(nlevel,ntime)
     253      real temp(nlevel,ntime),qv(nlevel,ntime),rh(nlevel,ntime)
     254      real theta(nlevel,ntime),rv(nlevel,ntime)
     255      real u(nlevel,ntime)
     256      real v(nlevel,ntime)
     257      real ug(nlevel,ntime)
     258      real vg(nlevel,ntime)
     259      real w(nlevel,ntime)
     260      real du(nlevel,ntime),hu(nlevel,ntime),vu(nlevel,ntime)
     261      real dv(nlevel,ntime),hv(nlevel,ntime),vv(nlevel,ntime)
     262      real dt(nlevel,ntime),ht(nlevel,ntime),vt(nlevel,ntime)
     263      real dtrad(nlevel,ntime)
     264      real dq(nlevel,ntime),hq(nlevel,ntime),vq(nlevel,ntime)
     265      real dth(nlevel,ntime),hth(nlevel,ntime),vth(nlevel,ntime)
     266      real dr(nlevel,ntime),hr(nlevel,ntime),vr(nlevel,ntime)
     267      real flat(ntime),sens(ntime),ts(ntime)
    194268
    195269
    196270      integer nid, ierr,rid
    197271      integer nbvar3d
    198       parameter(nbvar3d=30)
     272      parameter(nbvar3d=34)
    199273      integer var3didin(nbvar3d)
    200274
    … …  
    204278           stop 'lev'
    205279         endif
    206 
    207 
    208       ierr=NF_INQ_VARID(nid,"temp",var3didin(2))
     280     
     281      ierr=NF_INQ_VARID(nid,"pp",var3didin(2))
     282         if(ierr/=NF_NOERR) then
     283           write(*,*) NF_STRERROR(ierr)
     284           stop 'plev'
     285         endif
     286
     287
     288      ierr=NF_INQ_VARID(nid,"temp",var3didin(3))
    209289         if(ierr/=NF_NOERR) then
    210290           write(*,*) NF_STRERROR(ierr)
    … …  
    212292         endif
    213293
    214       ierr=NF_INQ_VARID(nid,"qv",var3didin(3))
     294      ierr=NF_INQ_VARID(nid,"qv",var3didin(4))
    215295         if(ierr/=NF_NOERR) then
    216296           write(*,*) NF_STRERROR(ierr)
    … …  
    218298         endif
    219299
    220       ierr=NF_INQ_VARID(nid,"u",var3didin(4))
     300      ierr=NF_INQ_VARID(nid,"rh",var3didin(5))
     301         if(ierr/=NF_NOERR) then
     302           write(*,*) NF_STRERROR(ierr)
     303           stop 'rh'
     304         endif
     305
     306      ierr=NF_INQ_VARID(nid,"theta",var3didin(6))
     307         if(ierr/=NF_NOERR) then
     308           write(*,*) NF_STRERROR(ierr)
     309           stop 'theta'
     310         endif
     311
     312      ierr=NF_INQ_VARID(nid,"rv",var3didin(7))
     313         if(ierr/=NF_NOERR) then
     314           write(*,*) NF_STRERROR(ierr)
     315           stop 'rv'
     316         endif
     317
     318
     319      ierr=NF_INQ_VARID(nid,"u",var3didin(8))
    221320         if(ierr/=NF_NOERR) then
    222321           write(*,*) NF_STRERROR(ierr)
    … …  
    224323         endif
    225324
    226       ierr=NF_INQ_VARID(nid,"v",var3didin(5))
     325      ierr=NF_INQ_VARID(nid,"v",var3didin(9))
    227326         if(ierr/=NF_NOERR) then
    228327           write(*,*) NF_STRERROR(ierr)
    … …  
    230329         endif
    231330
    232       ierr=NF_INQ_VARID(nid,"dw",var3didin(6))
    233          if(ierr/=NF_NOERR) then
    234            write(*,*) NF_STRERROR(ierr)
    235            stop 'dw'
    236          endif
    237 
    238       ierr=NF_INQ_VARID(nid,"dt",var3didin(7))
    239          if(ierr/=NF_NOERR) then
    240            write(*,*) NF_STRERROR(ierr)
    241            stop 'dt'
    242          endif
    243 
    244       ierr=NF_INQ_VARID(nid,"dq",var3didin(8))
    245          if(ierr/=NF_NOERR) then
    246            write(*,*) NF_STRERROR(ierr)
    247            stop 'dq'
     331       ierr=NF_INQ_VARID(nid,"ug",var3didin(10))
     332         if(ierr/=NF_NOERR) then
     333           write(*,*) NF_STRERROR(ierr)
     334           stop 'ug'
     335         endif
     336
     337      ierr=NF_INQ_VARID(nid,"vg",var3didin(11))
     338         if(ierr/=NF_NOERR) then
     339           write(*,*) NF_STRERROR(ierr)
     340           stop 'vg'
     341         endif
     342
     343      ierr=NF_INQ_VARID(nid,"w",var3didin(12))
     344         if(ierr/=NF_NOERR) then
     345           write(*,*) NF_STRERROR(ierr)
     346           stop 'w'
     347         endif
     348
     349      ierr=NF_INQ_VARID(nid,"advu",var3didin(13))
     350         if(ierr/=NF_NOERR) then
     351           write(*,*) NF_STRERROR(ierr)
     352           stop 'advu'
     353         endif
     354
     355      ierr=NF_INQ_VARID(nid,"hu",var3didin(14))
     356         if(ierr/=NF_NOERR) then
     357           write(*,*) NF_STRERROR(ierr)
     358           stop 'hu'
     359         endif
     360
     361       ierr=NF_INQ_VARID(nid,"vu",var3didin(15))
     362         if(ierr/=NF_NOERR) then
     363           write(*,*) NF_STRERROR(ierr)
     364           stop 'vu'
     365         endif
     366
     367       ierr=NF_INQ_VARID(nid,"advv",var3didin(16))
     368         if(ierr/=NF_NOERR) then
     369           write(*,*) NF_STRERROR(ierr)
     370           stop 'advv'
     371         endif
     372
     373      ierr=NF_INQ_VARID(nid,"hv",var3didin(17))
     374         if(ierr/=NF_NOERR) then
     375           write(*,*) NF_STRERROR(ierr)
     376           stop 'hv'
     377         endif
     378
     379       ierr=NF_INQ_VARID(nid,"vv",var3didin(18))
     380         if(ierr/=NF_NOERR) then
     381           write(*,*) NF_STRERROR(ierr)
     382           stop 'vv'
     383         endif
     384
     385      ierr=NF_INQ_VARID(nid,"advT",var3didin(19))
     386         if(ierr/=NF_NOERR) then
     387           write(*,*) NF_STRERROR(ierr)
     388           stop 'advT'
     389         endif
     390
     391      ierr=NF_INQ_VARID(nid,"hT",var3didin(20))
     392         if(ierr/=NF_NOERR) then
     393           write(*,*) NF_STRERROR(ierr)
     394           stop 'hT'
     395         endif
     396
     397      ierr=NF_INQ_VARID(nid,"vT",var3didin(21))
     398         if(ierr/=NF_NOERR) then
     399           write(*,*) NF_STRERROR(ierr)
     400           stop 'vT'
     401         endif
     402
     403      ierr=NF_INQ_VARID(nid,"advq",var3didin(22))
     404         if(ierr/=NF_NOERR) then
     405           write(*,*) NF_STRERROR(ierr)
     406           stop 'advq'
    248407         endif
    249408     
    250       ierr=NF_INQ_VARID(nid,"sens",var3didin(9))
     409      ierr=NF_INQ_VARID(nid,"hq",var3didin(23))
     410         if(ierr/=NF_NOERR) then
     411           write(*,*) NF_STRERROR(ierr)
     412           stop 'hq'
     413         endif
     414
     415      ierr=NF_INQ_VARID(nid,"vq",var3didin(24))
     416         if(ierr/=NF_NOERR) then
     417           write(*,*) NF_STRERROR(ierr)
     418           stop 'vq'
     419         endif
     420
     421      ierr=NF_INQ_VARID(nid,"advth",var3didin(25))
     422         if(ierr/=NF_NOERR) then
     423           write(*,*) NF_STRERROR(ierr)
     424           stop 'advth'
     425         endif
     426
     427      ierr=NF_INQ_VARID(nid,"hth",var3didin(26))
     428         if(ierr/=NF_NOERR) then
     429           write(*,*) NF_STRERROR(ierr)
     430           stop 'hth'
     431         endif
     432
     433      ierr=NF_INQ_VARID(nid,"vth",var3didin(27))
     434         if(ierr/=NF_NOERR) then
     435           write(*,*) NF_STRERROR(ierr)
     436           stop 'vth'
     437         endif
     438
     439      ierr=NF_INQ_VARID(nid,"advr",var3didin(28))
     440         if(ierr/=NF_NOERR) then
     441           write(*,*) NF_STRERROR(ierr)
     442           stop 'advr'
     443         endif
     444     
     445      ierr=NF_INQ_VARID(nid,"hr",var3didin(29))
     446         if(ierr/=NF_NOERR) then
     447           write(*,*) NF_STRERROR(ierr)
     448           stop 'hr'
     449         endif
     450
     451      ierr=NF_INQ_VARID(nid,"vr",var3didin(30))
     452         if(ierr/=NF_NOERR) then
     453           write(*,*) NF_STRERROR(ierr)
     454           stop 'vr'
     455         endif
     456
     457      ierr=NF_INQ_VARID(nid,"radT",var3didin(31))
     458         if(ierr/=NF_NOERR) then
     459           write(*,*) NF_STRERROR(ierr)
     460           stop 'radT'
     461         endif
     462
     463      ierr=NF_INQ_VARID(nid,"sens",var3didin(32))
    251464         if(ierr/=NF_NOERR) then
    252465           write(*,*) NF_STRERROR(ierr)
    … …  
    254467         endif
    255468
    256       ierr=NF_INQ_VARID(nid,"flat",var3didin(10))
     469      ierr=NF_INQ_VARID(nid,"flat",var3didin(33))
    257470         if(ierr/=NF_NOERR) then
    258471           write(*,*) NF_STRERROR(ierr)
    … …  
    260473         endif
    261474
    262       ierr=NF_INQ_VARID(nid,"pp",var3didin(11))
    263          if(ierr/=NF_NOERR) then
    264            write(*,*) NF_STRERROR(ierr)
    265       endif
    266 
    267 !dimensions lecture
    268 !      call catchaxis(nid,ntime,nlevel,time,z,ierr)
     475      ierr=NF_INQ_VARID(nid,"ts",var3didin(34))
     476         if(ierr/=NF_NOERR) then
     477           write(*,*) NF_STRERROR(ierr)
     478           stop 'ts'
     479         endif
    269480 
    270481#ifdef NC_DOUBLE
    … …  
    280491
    281492#ifdef NC_DOUBLE
    282          ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),temp)
    283 #else
    284          ierr = NF_GET_VAR_REAL(nid,var3didin(2),temp)
    285 #endif
    286          if(ierr/=NF_NOERR) then
    287             write(*,*) NF_STRERROR(ierr)
    288             stop "getvarup"
    289          endif
    290 !          write(*,*)'lecture th ok',temp
    291 
    292 #ifdef NC_DOUBLE
    293          ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qv)
    294 #else
    295          ierr = NF_GET_VAR_REAL(nid,var3didin(3),qv)
     493         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),pp)
     494#else
     495         ierr = NF_GET_VAR_REAL(nid,var3didin(2),pp)
     496#endif
     497         if(ierr/=NF_NOERR) then
     498            write(*,*) NF_STRERROR(ierr)
     499            stop "getvarup"
     500         endif
     501!          write(*,*)'lecture pp ok',pp
     502
     503
     504#ifdef NC_DOUBLE
     505         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),temp)
     506#else
     507         ierr = NF_GET_VAR_REAL(nid,var3didin(3),temp)
     508#endif
     509         if(ierr/=NF_NOERR) then
     510            write(*,*) NF_STRERROR(ierr)
     511            stop "getvarup"
     512         endif
     513!          write(*,*)'lecture T ok',temp
     514
     515#ifdef NC_DOUBLE
     516         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),qv)
     517#else
     518         ierr = NF_GET_VAR_REAL(nid,var3didin(4),qv)
    296519#endif
    297520         if(ierr/=NF_NOERR) then
    … …  
    302525 
    303526#ifdef NC_DOUBLE
    304          ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)
    305 #else
    306          ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)
     527         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),rh)
     528#else
     529         ierr = NF_GET_VAR_REAL(nid,var3didin(5),rh)
     530#endif
     531         if(ierr/=NF_NOERR) then
     532            write(*,*) NF_STRERROR(ierr)
     533            stop "getvarup"
     534         endif
     535!          write(*,*)'lecture rh ok',rh
     536
     537#ifdef NC_DOUBLE
     538         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),theta)
     539#else
     540         ierr = NF_GET_VAR_REAL(nid,var3didin(6),theta)
     541#endif
     542         if(ierr/=NF_NOERR) then
     543            write(*,*) NF_STRERROR(ierr)
     544            stop "getvarup"
     545         endif
     546!          write(*,*)'lecture theta ok',theta
     547
     548#ifdef NC_DOUBLE
     549         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),rv)
     550#else
     551         ierr = NF_GET_VAR_REAL(nid,var3didin(7),rv)
     552#endif
     553         if(ierr/=NF_NOERR) then
     554            write(*,*) NF_STRERROR(ierr)
     555            stop "getvarup"
     556         endif
     557!          write(*,*)'lecture rv ok',rv
     558
     559#ifdef NC_DOUBLE
     560         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),u)
     561#else
     562         ierr = NF_GET_VAR_REAL(nid,var3didin(8),u)
    307563#endif
    308564         if(ierr/=NF_NOERR) then
    … …  
    313569
    314570#ifdef NC_DOUBLE
    315          ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)
    316 #else
    317          ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)
     571         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),v)
     572#else
     573         ierr = NF_GET_VAR_REAL(nid,var3didin(9),v)
    318574#endif
    319575         if(ierr/=NF_NOERR) then
    … …  
    324580
    325581#ifdef NC_DOUBLE
    326          ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),dw)
    327 #else
    328          ierr = NF_GET_VAR_REAL(nid,var3didin(6),dw)
    329 #endif
    330          if(ierr/=NF_NOERR) then
    331             write(*,*) NF_STRERROR(ierr)
    332             stop "getvarup"
    333          endif
    334 !          write(*,*)'lecture w ok',dw
    335 
    336 #ifdef NC_DOUBLE
    337          ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),dt)
    338 #else
    339          ierr = NF_GET_VAR_REAL(nid,var3didin(7),dt)
     582         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),ug)
     583#else
     584         ierr = NF_GET_VAR_REAL(nid,var3didin(10),ug)
     585#endif
     586         if(ierr/=NF_NOERR) then
     587            write(*,*) NF_STRERROR(ierr)
     588            stop "getvarup"
     589         endif
     590!          write(*,*)'lecture ug ok',ug
     591
     592#ifdef NC_DOUBLE
     593         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),vg)
     594#else
     595         ierr = NF_GET_VAR_REAL(nid,var3didin(11),vg)
     596#endif
     597         if(ierr/=NF_NOERR) then
     598            write(*,*) NF_STRERROR(ierr)
     599            stop "getvarup"
     600         endif
     601!          write(*,*)'lecture vg ok',vg
     602
     603#ifdef NC_DOUBLE
     604         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),w)
     605#else
     606         ierr = NF_GET_VAR_REAL(nid,var3didin(12),w)
     607#endif
     608         if(ierr/=NF_NOERR) then
     609            write(*,*) NF_STRERROR(ierr)
     610            stop "getvarup"
     611         endif
     612!          write(*,*)'lecture w ok',w
     613
     614#ifdef NC_DOUBLE
     615         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),du)
     616#else
     617         ierr = NF_GET_VAR_REAL(nid,var3didin(13),du)
     618#endif
     619         if(ierr/=NF_NOERR) then
     620            write(*,*) NF_STRERROR(ierr)
     621            stop "getvarup"
     622         endif
     623!          write(*,*)'lecture du ok',du
     624
     625#ifdef NC_DOUBLE
     626         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),hu)
     627#else
     628         ierr = NF_GET_VAR_REAL(nid,var3didin(14),hu)
     629#endif
     630         if(ierr/=NF_NOERR) then
     631            write(*,*) NF_STRERROR(ierr)
     632            stop "getvarup"
     633         endif
     634!          write(*,*)'lecture hu ok',hu
     635
     636#ifdef NC_DOUBLE
     637         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),vu)
     638#else
     639         ierr = NF_GET_VAR_REAL(nid,var3didin(15),vu)
     640#endif
     641         if(ierr/=NF_NOERR) then
     642            write(*,*) NF_STRERROR(ierr)
     643            stop "getvarup"
     644         endif
     645!          write(*,*)'lecture vu ok',vu
     646
     647#ifdef NC_DOUBLE
     648         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),dv)
     649#else
     650         ierr = NF_GET_VAR_REAL(nid,var3didin(16),dv)
     651#endif
     652         if(ierr/=NF_NOERR) then
     653            write(*,*) NF_STRERROR(ierr)
     654            stop "getvarup"
     655         endif
     656!          write(*,*)'lecture dv ok',dv
     657
     658#ifdef NC_DOUBLE
     659         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),hv)
     660#else
     661         ierr = NF_GET_VAR_REAL(nid,var3didin(17),hv)
     662#endif
     663         if(ierr/=NF_NOERR) then
     664            write(*,*) NF_STRERROR(ierr)
     665            stop "getvarup"
     666         endif
     667!          write(*,*)'lecture hv ok',hv
     668
     669#ifdef NC_DOUBLE
     670         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),vv)
     671#else
     672         ierr = NF_GET_VAR_REAL(nid,var3didin(18),vv)
     673#endif
     674         if(ierr/=NF_NOERR) then
     675            write(*,*) NF_STRERROR(ierr)
     676            stop "getvarup"
     677         endif
     678!          write(*,*)'lecture vv ok',vv
     679
     680#ifdef NC_DOUBLE
     681         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),dt)
     682#else
     683         ierr = NF_GET_VAR_REAL(nid,var3didin(19),dt)
    340684#endif
    341685         if(ierr/=NF_NOERR) then
    … …  
    346690
    347691#ifdef NC_DOUBLE
    348          ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),dq)
    349 #else
    350          ierr = NF_GET_VAR_REAL(nid,var3didin(8),dq)
     692         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),ht)
     693#else
     694         ierr = NF_GET_VAR_REAL(nid,var3didin(20),ht)
     695#endif
     696         if(ierr/=NF_NOERR) then
     697            write(*,*) NF_STRERROR(ierr)
     698            stop "getvarup"
     699         endif
     700!          write(*,*)'lecture ht ok',ht
     701
     702#ifdef NC_DOUBLE
     703         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(21),vt)
     704#else
     705         ierr = NF_GET_VAR_REAL(nid,var3didin(21),vt)
     706#endif
     707         if(ierr/=NF_NOERR) then
     708            write(*,*) NF_STRERROR(ierr)
     709            stop "getvarup"
     710         endif
     711!          write(*,*)'lecture vt ok',vt
     712
     713#ifdef NC_DOUBLE
     714         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(22),dq)
     715#else
     716         ierr = NF_GET_VAR_REAL(nid,var3didin(22),dq)
    351717#endif
    352718         if(ierr/=NF_NOERR) then
    … …  
    357723
    358724#ifdef NC_DOUBLE
    359          ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),sens)
    360 #else
    361          ierr = NF_GET_VAR_REAL(nid,var3didin(9),sens)
     725         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(23),hq)
     726#else
     727         ierr = NF_GET_VAR_REAL(nid,var3didin(23),hq)
     728#endif
     729         if(ierr/=NF_NOERR) then
     730            write(*,*) NF_STRERROR(ierr)
     731            stop "getvarup"
     732         endif
     733!          write(*,*)'lecture hq ok',hq
     734
     735#ifdef NC_DOUBLE
     736         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(24),vq)
     737#else
     738         ierr = NF_GET_VAR_REAL(nid,var3didin(24),vq)
     739#endif
     740         if(ierr/=NF_NOERR) then
     741            write(*,*) NF_STRERROR(ierr)
     742            stop "getvarup"
     743         endif
     744!          write(*,*)'lecture vq ok',vq
     745
     746#ifdef NC_DOUBLE
     747         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(25),dth)
     748#else
     749         ierr = NF_GET_VAR_REAL(nid,var3didin(25),dth)
     750#endif
     751         if(ierr/=NF_NOERR) then
     752            write(*,*) NF_STRERROR(ierr)
     753            stop "getvarup"
     754         endif
     755!          write(*,*)'lecture dth ok',dth
     756
     757#ifdef NC_DOUBLE
     758         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(26),hth)
     759#else
     760         ierr = NF_GET_VAR_REAL(nid,var3didin(26),hth)
     761#endif
     762         if(ierr/=NF_NOERR) then
     763            write(*,*) NF_STRERROR(ierr)
     764            stop "getvarup"
     765         endif
     766!          write(*,*)'lecture hth ok',hth
     767
     768#ifdef NC_DOUBLE
     769         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(27),vth)
     770#else
     771         ierr = NF_GET_VAR_REAL(nid,var3didin(27),vth)
     772#endif
     773         if(ierr/=NF_NOERR) then
     774            write(*,*) NF_STRERROR(ierr)
     775            stop "getvarup"
     776         endif
     777!          write(*,*)'lecture vth ok',vth
     778
     779#ifdef NC_DOUBLE
     780         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(28),dr)
     781#else
     782         ierr = NF_GET_VAR_REAL(nid,var3didin(28),dr)
     783#endif
     784         if(ierr/=NF_NOERR) then
     785            write(*,*) NF_STRERROR(ierr)
     786            stop "getvarup"
     787         endif
     788!          write(*,*)'lecture dr ok',dr
     789
     790#ifdef NC_DOUBLE
     791         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(29),hr)
     792#else
     793         ierr = NF_GET_VAR_REAL(nid,var3didin(29),hr)
     794#endif
     795         if(ierr/=NF_NOERR) then
     796            write(*,*) NF_STRERROR(ierr)
     797            stop "getvarup"
     798         endif
     799!          write(*,*)'lecture hr ok',hr
     800
     801#ifdef NC_DOUBLE
     802         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(30),vr)
     803#else
     804         ierr = NF_GET_VAR_REAL(nid,var3didin(30),vr)
     805#endif
     806         if(ierr/=NF_NOERR) then
     807            write(*,*) NF_STRERROR(ierr)
     808            stop "getvarup"
     809         endif
     810!          write(*,*)'lecture vr ok',vr
     811
     812#ifdef NC_DOUBLE
     813         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(31),dtrad)
     814#else
     815         ierr = NF_GET_VAR_REAL(nid,var3didin(31),dtrad)
     816#endif
     817         if(ierr/=NF_NOERR) then
     818            write(*,*) NF_STRERROR(ierr)
     819            stop "getvarup"
     820         endif
     821!          write(*,*)'lecture dtrad ok',dtrad
     822
     823#ifdef NC_DOUBLE
     824         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(32),sens)
     825#else
     826         ierr = NF_GET_VAR_REAL(nid,var3didin(32),sens)
    362827#endif
    363828         if(ierr/=NF_NOERR) then
    … …  
    368833
    369834#ifdef NC_DOUBLE
    370          ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),flat)
    371 #else
    372          ierr = NF_GET_VAR_REAL(nid,var3didin(10),flat)
     835         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(33),flat)
     836#else
     837         ierr = NF_GET_VAR_REAL(nid,var3didin(33),flat)
    373838#endif
    374839         if(ierr/=NF_NOERR) then
    … …  
    379844
    380845#ifdef NC_DOUBLE
    381          ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pp)
    382 #else
    383          ierr = NF_GET_VAR_REAL(nid,var3didin(11),pp)
    384 #endif
    385          if(ierr/=NF_NOERR) then
    386             write(*,*) NF_STRERROR(ierr)
    387             stop "getvarup"
    388          endif
    389 !          write(*,*)'lecture pp ok',pp
     846         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(34),ts)
     847#else
     848         ierr = NF_GET_VAR_REAL(nid,var3didin(34),ts)
     849#endif
     850         if(ierr/=NF_NOERR) then
     851            write(*,*) NF_STRERROR(ierr)
     852            stop "getvarup"
     853         endif
     854!          write(*,*)'lecture ts ok',ts
    390855
    391856         return
    392          end subroutine read_amma
     857         end subroutine read_cas
    393858!======================================================================
    394         SUBROUTINE interp_amma_time(day,day1,annee_ref                     &
    395      &         ,year_ini_amma,day_ini_amma,nt_amma,dt_amma,nlev_amma       &
    396      &         ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma               &
    397      &         ,vitw_prof,ht_prof,hq_prof,lat_prof,sens_prof)
     859        SUBROUTINE interp_case_time(day,day1,annee_ref                &
     860     &         ,year_ini_cas,day_ini_cas,nt_cas,pdt_cas,nlev_cas      &
     861     &         ,ts_cas,plev_cas,t_cas,q_cas,u_cas,v_cas               &
     862     &         ,ug_cas,vg_cas,vitw_cas,du_cas,hu_cas,vu_cas           &
     863     &         ,dv_cas,hv_cas,vv_cas,dt_cas,ht_cas,vt_cas,dtrad_cas   &
     864     &         ,dq_cas,hq_cas,vq_cas,lat_cas,sens_cas          &
     865     &         ,ts_prof_cas,plev_prof_cas,t_prof_cas,q_prof_cas       &
     866     &         ,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas         &
     867     &         ,vitw_prof_cas,du_prof_cas,hu_prof_cas,vu_prof_cas     &
     868     &         ,dv_prof_cas,hv_prof_cas,vv_prof_cas,dt_prof_cas       &
     869     &         ,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas    &
     870     &         ,hq_prof_cas,vq_prof_cas,lat_prof_cas,sens_prof_cas)
     871         
     872
    398873        implicit none
    399874
    … …  
    403878! day: current julian day (e.g. 717538.2)
    404879! day1: first day of the simulation
    405 ! nt_amma: total nb of data in the forcing (e.g. 48 for AMMA)
    406 ! dt_amma: total time interval (in sec) between 2 forcing data (e.g. 30min for AMMA)
     880! nt_cas: total nb of data in the forcing
     881! pdt_cas: total time interval (in sec) between 2 forcing data
    407882!---------------------------------------------------------------------------------------
    408883
    … …  
    411886! inputs:
    412887        integer annee_ref
    413         integer nt_amma,nlev_amma
    414         integer year_ini_amma
    415         real day, day1,day_ini_amma,dt_amma
    416         real vitw_amma(nlev_amma,nt_amma)
    417         real ht_amma(nlev_amma,nt_amma)
    418         real hq_amma(nlev_amma,nt_amma)
    419         real lat_amma(nt_amma)
    420         real sens_amma(nt_amma)
     888        integer nt_cas,nlev_cas
     889        integer year_ini_cas
     890        real day_ini_cas
     891        real day, day1,pdt_cas
     892        real ts_cas(nt_cas)
     893        real plev_cas(nlev_cas,nt_cas)
     894        real t_cas(nlev_cas,nt_cas),q_cas(nlev_cas,nt_cas)
     895        real u_cas(nlev_cas,nt_cas),v_cas(nlev_cas,nt_cas)
     896        real ug_cas(nlev_cas,nt_cas),vg_cas(nlev_cas,nt_cas)
     897        real vitw_cas(nlev_cas,nt_cas)
     898        real du_cas(nlev_cas,nt_cas),hu_cas(nlev_cas,nt_cas),vu_cas(nlev_cas,nt_cas)
     899        real dv_cas(nlev_cas,nt_cas),hv_cas(nlev_cas,nt_cas),vv_cas(nlev_cas,nt_cas)
     900        real dt_cas(nlev_cas,nt_cas),ht_cas(nlev_cas,nt_cas),vt_cas(nlev_cas,nt_cas)
     901        real dtrad_cas(nlev_cas,nt_cas)
     902        real dq_cas(nlev_cas,nt_cas),hq_cas(nlev_cas,nt_cas),vq_cas(nlev_cas,nt_cas)
     903        real lat_cas(nt_cas)
     904        real sens_cas(nt_cas)
     905
    421906! outputs:
    422         real vitw_prof(nlev_amma)
    423         real ht_prof(nlev_amma)
    424         real hq_prof(nlev_amma)
    425         real lat_prof,sens_prof
     907        real plev_prof_cas(nlev_cas)
     908        real t_prof_cas(nlev_cas),q_prof_cas(nlev_cas)
     909        real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)
     910        real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas)
     911        real vitw_prof_cas(nlev_cas)
     912        real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)
     913        real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)
     914        real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas)
     915        real dtrad_prof_cas(nlev_cas)
     916        real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)
     917        real lat_prof_cas,sens_prof_cas,ts_prof_cas
    426918! local:
    427         integer it_amma1, it_amma2,k
    428         real timeit,time_amma1,time_amma2,frac
    429 
    430 
    431         if (forcing_type.eq.6) then
     919        integer it_cas1, it_cas2,k
     920        real timeit,time_cas1,time_cas2,frac
     921
     922
     923        print*,'Check time',day1,day_ini_cas,day_ini_cas+90
     924
     925        if ((forcing_type.eq.10).and.(1.eq.1)) then
     926! Check that initial day of the simulation consistent with the case:
     927       if (annee_ref.ne.2011) then
     928        print*,'Pour CINDY, annee_ref doit etre 2011'
     929        print*,'Changer annee_ref dans run.def'
     930        stop
     931       endif
     932       if (annee_ref.eq.2011 .and. day1.lt.day_ini_cas) then
     933        print*,'CINDY a débuté le 1 octobre 2011',day1,day_ini_cas
     934        print*,'Changer dayref dans run.def'
     935        stop
     936       endif
     937       if (annee_ref.eq.2011 .and. day1.gt.day_ini_cas+90) then
     938        print*,'CINDY a fini le 31 decembre'
     939        print*,'Changer dayref ou nday dans run.def',day1,day_ini_cas+90
     940        stop
     941       endif
     942       endif
     943!CR:ATTENTION TEST AMMA
     944      if ((forcing_type.eq.10).and.(1.eq.0)) then
    432945! Check that initial day of the simulation consistent with AMMA case:
    433946       if (annee_ref.ne.2006) then
    … …  
    436949        stop
    437950       endif
    438        if (annee_ref.eq.2006 .and. day1.lt.day_ini_amma) then
    439         print*,'AMMA a débuté le 10 juillet 2006',day1,day_ini_amma
     951       if (annee_ref.eq.2006 .and. day1.lt.day_ini_cas) then
     952        print*,'AMMA a débuté le 10 juillet 2006',day1,day_ini_cas
    440953        print*,'Changer dayref dans run.def'
    441954        stop
    442955       endif
    443        if (annee_ref.eq.2006 .and. day1.gt.day_ini_amma+1) then
     956       if (annee_ref.eq.2006 .and. day1.gt.day_ini_cas+1) then
    444957        print*,'AMMA a fini le 11 juillet'
    445958        print*,'Changer dayref ou nday dans run.def'
    … …  
    448961       endif
    449962
    450 ! Determine timestep relative to the 1st day of AMMA:
     963! Determine timestep relative to the 1st day:
    451964!       timeit=(day-day1)*86400.
    452965!       if (annee_ref.eq.1992) then
    453 !        timeit=(day-day_ini_toga)*86400.
     966!        timeit=(day-day_ini_cas)*86400.
    454967!       else
    455968!        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
    456969!       endif
    457       timeit=(day-day_ini_amma)*86400
     970      timeit=(day-day_ini_cas)*86400
    458971
    459972! Determine the closest observation times:
    460 !       it_amma1=INT(timeit/dt_amma)+1
    461 !       it_amma2=it_amma1 + 1
    462 !       time_amma1=(it_amma1-1)*dt_amma
    463 !       time_amma2=(it_amma2-1)*dt_amma
    464 
    465        it_amma1=INT(timeit/dt_amma)+1
    466        IF (it_amma1 .EQ. nt_amma) THEN
    467        it_amma2=it_amma1
     973!       it_cas1=INT(timeit/pdt_cas)+1
     974!       it_cas2=it_cas1 + 1
     975!       time_cas1=(it_cas1-1)*pdt_cas
     976!       time_cas2=(it_cas2-1)*pdt_cas
     977
     978       it_cas1=INT(timeit/pdt_cas)+1
     979       IF (it_cas1 .EQ. nt_cas) THEN
     980       it_cas2=it_cas1
    468981       ELSE
    469        it_amma2=it_amma1 + 1
     982       it_cas2=it_cas1 + 1
    470983       ENDIF
    471        time_amma1=(it_amma1-1)*dt_amma
    472        time_amma2=(it_amma2-1)*dt_amma
    473 
    474        if (it_amma1 .gt. nt_amma) then
    475         write(*,*) 'PB-stop: day, it_amma1, it_amma2, timeit: '            &
    476      &        ,day,day_ini_amma,it_amma1,it_amma2,timeit/86400.
     984       time_cas1=(it_cas1-1)*pdt_cas
     985       time_cas2=(it_cas2-1)*pdt_cas
     986
     987       if (it_cas1 .gt. nt_cas) then
     988        write(*,*) 'PB-stop: day, it_cas1, it_cas2, timeit: '            &
     989     &        ,day,day_ini_cas,it_cas1,it_cas2,timeit/86400.
    477990        stop
    478991       endif
    479992
    480993! time interpolation:
    481        frac=(time_amma2-timeit)/(time_amma2-time_amma1)
     994       frac=(time_cas2-timeit)/(time_cas2-time_cas1)
    482995       frac=max(frac,0.0)
    483996
    484        lat_prof = lat_amma(it_amma2)                                       &
    485      &          -frac*(lat_amma(it_amma2)-lat_amma(it_amma1))
    486        sens_prof = sens_amma(it_amma2)                                     &
    487      &          -frac*(sens_amma(it_amma2)-sens_amma(it_amma1))
    488 
    489        do k=1,nlev_amma
    490         vitw_prof(k) = vitw_amma(k,it_amma2)                               &
    491      &          -frac*(vitw_amma(k,it_amma2)-vitw_amma(k,it_amma1))
    492         ht_prof(k) = ht_amma(k,it_amma2)                                   &
    493      &          -frac*(ht_amma(k,it_amma2)-ht_amma(k,it_amma1))
    494         hq_prof(k) = hq_amma(k,it_amma2)                                   &
    495      &          -frac*(hq_amma(k,it_amma2)-hq_amma(k,it_amma1))
     997       lat_prof_cas = lat_cas(it_cas2)                                       &
     998     &          -frac*(lat_cas(it_cas2)-lat_cas(it_cas1))
     999       sens_prof_cas = sens_cas(it_cas2)                                     &
     1000     &          -frac*(sens_cas(it_cas2)-sens_cas(it_cas1))
     1001       ts_prof_cas = ts_cas(it_cas2)                                     &
     1002     &          -frac*(ts_cas(it_cas2)-ts_cas(it_cas1))
     1003
     1004       do k=1,nlev_cas
     1005        plev_prof_cas(k) = plev_cas(k,it_cas2)                               &
     1006     &          -frac*(plev_cas(k,it_cas2)-plev_cas(k,it_cas1))
     1007        t_prof_cas(k) = t_cas(k,it_cas2)                               &
     1008     &          -frac*(t_cas(k,it_cas2)-t_cas(k,it_cas1))
     1009        q_prof_cas(k) = q_cas(k,it_cas2)                               &
     1010     &          -frac*(q_cas(k,it_cas2)-q_cas(k,it_cas1))
     1011        u_prof_cas(k) = u_cas(k,it_cas2)                               &
     1012     &          -frac*(u_cas(k,it_cas2)-u_cas(k,it_cas1))
     1013        v_prof_cas(k) = v_cas(k,it_cas2)                               &
     1014     &          -frac*(v_cas(k,it_cas2)-v_cas(k,it_cas1))
     1015        ug_prof_cas(k) = ug_cas(k,it_cas2)                               &
     1016     &          -frac*(ug_cas(k,it_cas2)-ug_cas(k,it_cas1))
     1017        vg_prof_cas(k) = vg_cas(k,it_cas2)                               &
     1018     &          -frac*(vg_cas(k,it_cas2)-vg_cas(k,it_cas1))
     1019        vitw_prof_cas(k) = vitw_cas(k,it_cas2)                               &
     1020     &          -frac*(vitw_cas(k,it_cas2)-vitw_cas(k,it_cas1))
     1021        du_prof_cas(k) = du_cas(k,it_cas2)                                   &
     1022     &          -frac*(du_cas(k,it_cas2)-du_cas(k,it_cas1))
     1023        hu_prof_cas(k) = hu_cas(k,it_cas2)                                   &
     1024     &          -frac*(hu_cas(k,it_cas2)-hu_cas(k,it_cas1))
     1025        vu_prof_cas(k) = vu_cas(k,it_cas2)                                   &
     1026     &          -frac*(vu_cas(k,it_cas2)-vu_cas(k,it_cas1))
     1027        dv_prof_cas(k) = dv_cas(k,it_cas2)                                   &
     1028     &          -frac*(dv_cas(k,it_cas2)-dv_cas(k,it_cas1))
     1029        hv_prof_cas(k) = hv_cas(k,it_cas2)                                   &
     1030     &          -frac*(hv_cas(k,it_cas2)-hv_cas(k,it_cas1))
     1031        vv_prof_cas(k) = vv_cas(k,it_cas2)                                   &
     1032     &          -frac*(vv_cas(k,it_cas2)-vv_cas(k,it_cas1))
     1033        dt_prof_cas(k) = dt_cas(k,it_cas2)                                   &
     1034     &          -frac*(dt_cas(k,it_cas2)-dt_cas(k,it_cas1))
     1035        ht_prof_cas(k) = ht_cas(k,it_cas2)                                   &
     1036     &          -frac*(ht_cas(k,it_cas2)-ht_cas(k,it_cas1))
     1037        vt_prof_cas(k) = vt_cas(k,it_cas2)                                   &
     1038     &          -frac*(vt_cas(k,it_cas2)-vt_cas(k,it_cas1))
     1039        dtrad_prof_cas(k) = dtrad_cas(k,it_cas2)                                   &
     1040     &          -frac*(dtrad_cas(k,it_cas2)-dtrad_cas(k,it_cas1))
     1041        dq_prof_cas(k) = dq_cas(k,it_cas2)                                   &
     1042     &          -frac*(dq_cas(k,it_cas2)-dq_cas(k,it_cas1))
     1043        hq_prof_cas(k) = hq_cas(k,it_cas2)                                   &
     1044     &          -frac*(hq_cas(k,it_cas2)-hq_cas(k,it_cas1))
     1045        vq_prof_cas(k) = vq_cas(k,it_cas2)                                   &
     1046     &          -frac*(vq_cas(k,it_cas2)-vq_cas(k,it_cas1))
    4961047        enddo
    4971048
    … …  
    4991050        END
    5001051
     1052!**********************************************************************************************

  • LMDZ5/branches/testing/libf/phylmd/nuage.h

    r2056 r2220  
    55      REAL exposant_glace
    66      REAL rei_min,rei_max
     7      REAL tau_cld_cv,coefw_cld_cv
    78
    8       INTEGER iflag_t_glace
     9      INTEGER iflag_t_glace,iflag_cld_cv
    910
    1011      common /nuagecom/ rad_froid,rad_chau1, rad_chau2,t_glace_max,     &
    1112     &                  t_glace_min,exposant_glace,rei_min,rei_max,     &
    12      &                  iflag_t_glace
     13     &                  tau_cld_cv,coefw_cld_cv,                        &
     14     &                  iflag_t_glace,iflag_cld_cv
    1315!$OMP THREADPRIVATE(/nuagecom/)

  • LMDZ5/branches/testing/libf/phylmd/ocean_slab_mod.F90

    r2073 r2220  
    88  USE dimphy
    99  USE indice_sol_mod
     10  USE surface_data
    1011
    1112  IMPLICIT NONE
    1213  PRIVATE
    13   PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice!, ocean_slab_ice
     14  PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice
     15
     16!****************************************************************************************
     17! Global saved variables
     18!****************************************************************************************
    1419
    1520  INTEGER, PRIVATE, SAVE                           :: cpl_pas
    … …  
    2126  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: tslab
    2227  !$OMP THREADPRIVATE(tslab)
    23   REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE  :: pctsrf
    24   !$OMP THREADPRIVATE(pctsrf)
     28  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: fsic
     29  !$OMP THREADPRIVATE(fsic)
     30  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: tice
     31  !$OMP THREADPRIVATE(tice)
     32  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: seaice
     33  !$OMP THREADPRIVATE(seaice)
    2534  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bils
    2635  !$OMP THREADPRIVATE(slab_bils)
    2736  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bils_cum
    2837  !$OMP THREADPRIVATE(bils_cum)
     38  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bilg
     39  !$OMP THREADPRIVATE(slab_bilg)
     40  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bilg_cum
     41  !$OMP THREADPRIVATE(bilg_cum)
     42
     43!****************************************************************************************
     44! Parameters (could be read in def file: move to slab_init)
     45!****************************************************************************************
     46! snow and ice physical characteristics:
     47    REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp
     48    REAL, PARAMETER :: t_melt=273.15   ! melting ice temp
     49    REAL, PARAMETER :: sno_den=300. !mean snow density, kg/m3
     50    REAL, PARAMETER :: ice_den=917.
     51    REAL, PARAMETER :: sea_den=1025.
     52    REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity
     53    REAL, PARAMETER :: sno_cond=0.31*sno_den
     54    REAL, PARAMETER :: ice_cap=2067.   ! specific heat capacity, snow and ice
     55    REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice
     56
     57! control of snow and ice cover & freeze / melt (heights converted to kg/m2)
     58    REAL, PARAMETER :: snow_min=0.05*sno_den !critical snow height 5 cm
     59    REAL, PARAMETER :: snow_wfact=0.4 ! max fraction of falling snow blown into ocean
     60    REAL, PARAMETER :: ice_frac_min=0.001
     61    REAL, PARAMETER :: ice_frac_max=1. ! less than 1. if min leads fraction
     62    REAL, PARAMETER :: h_ice_min=0.01*ice_den ! min ice thickness
     63    REAL, PARAMETER :: h_ice_thin=0.15*ice_den ! thin ice thickness
     64    ! below ice_thin, priority is melt lateral / grow height
     65    ! ice_thin is also height of new ice
     66    REAL, PARAMETER :: h_ice_thick=2.5*ice_den ! thin ice thickness
     67    ! above ice_thick, priority is melt height / grow lateral
     68    REAL, PARAMETER :: h_ice_new=1.*ice_den ! max height of new open ocean ice
     69    REAL, PARAMETER :: h_ice_max=10.*ice_den ! max ice height
     70
     71! albedo  and radiation parameters
     72    REAL, PARAMETER :: alb_sno_min=0.55 !min snow albedo
     73    REAL, PARAMETER :: alb_sno_del=0.3  !max snow albedo = min + del
     74    REAL, PARAMETER :: alb_ice_dry=0.75 !dry thick ice
     75    REAL, PARAMETER :: alb_ice_wet=0.66 !melting thick ice
     76    REAL, PARAMETER :: pen_frac=0.3 !fraction of penetrating shortwave (no snow)
     77    REAL, PARAMETER :: pen_ext=1.5 !extinction of penetrating shortwave (m-1)
     78
     79!****************************************************************************************
    2980
    3081CONTAINS
    … …  
    56107! Allocate surface fraction read from restart file
    57108!****************************************************************************************
    58     ALLOCATE(pctsrf(klon,nbsrf), stat = error)
     109    ALLOCATE(fsic(klon), stat = error)
    59110    IF (error /= 0) THEN
    60111       abort_message='Pb allocation tmp_pctsrf_slab'
    61112       CALL abort_gcm(modname,abort_message,1)
    62113    ENDIF
    63     pctsrf(:,:) = pctsrf_rst(:,:)
    64 
    65 !****************************************************************************************
    66 ! Allocate local variables
    67 !****************************************************************************************
     114    fsic(:)=0.
     115    WHERE (1.-zmasq(:)>EPSFRA)
     116        fsic(:) = pctsrf_rst(:,is_sic)/(1.-zmasq(:))
     117    END WHERE
     118
     119!****************************************************************************************
     120! Allocate saved variables
     121!****************************************************************************************
     122    ALLOCATE(tslab(klon,nslay), stat=error)
     123       IF (error /= 0) CALL abort_gcm &
     124         (modname,'pb allocation tslab', 1)
     125
    68126    ALLOCATE(slab_bils(klon), stat = error)
    69127    IF (error /= 0) THEN
    … …  
    79137    bils_cum(:) = 0.0   
    80138
     139    IF (version_ocean=='sicINT') THEN
     140        ALLOCATE(slab_bilg(klon), stat = error)
     141        IF (error /= 0) THEN
     142           abort_message='Pb allocation slab_bilg'
     143           CALL abort_gcm(modname,abort_message,1)
     144        ENDIF
     145        slab_bilg(:) = 0.0   
     146        ALLOCATE(bilg_cum(klon), stat = error)
     147        IF (error /= 0) THEN
     148           abort_message='Pb allocation slab_bilg_cum'
     149           CALL abort_gcm(modname,abort_message,1)
     150        ENDIF
     151        bilg_cum(:) = 0.0   
     152        ALLOCATE(tice(klon), stat = error)
     153        IF (error /= 0) THEN
     154           abort_message='Pb allocation slab_tice'
     155           CALL abort_gcm(modname,abort_message,1)
     156        ENDIF
     157        ALLOCATE(seaice(klon), stat = error)
     158        IF (error /= 0) THEN
     159           abort_message='Pb allocation slab_seaice'
     160           CALL abort_gcm(modname,abort_message,1)
     161        ENDIF
     162    END IF
     163
     164!****************************************************************************************
     165! Define some parameters
     166!****************************************************************************************
    81167! Layer thickness
    82168    ALLOCATE(slabh(nslay), stat = error)
    … …  
    94180    CALL getin('cpl_pas',cpl_pas)
    95181    print *,'cpl_pas',cpl_pas
     182
    96183 END SUBROUTINE ocean_slab_init
    97184!
    … …  
    101188
    102189    USE limit_read_mod
    103     USE surface_data
    104190
    105191!    INCLUDE "clesphys.h"
    … …  
    119205       CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, is_modified)
    120206    ELSE
    121        pctsrf_chg(:,:)=pctsrf(:,:)
     207       pctsrf_chg(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:))
     208       pctsrf_chg(:,is_sic)=fsic(:)*(1.-zmasq(:))
    122209       is_modified=.TRUE.
    123210    END IF
    … …  
    133220       AcoefU, AcoefV, BcoefU, BcoefV, &
    134221       ps, u1, v1, tsurf_in, &
    135        radsol, snow, agesno, &
     222       radsol, snow, &
    136223       qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
    137224       tsurf_new, dflux_s, dflux_l, qflux)
    138225   
    139226    USE calcul_fluxs_mod
    140     USE surface_data
    141227
    142228    INCLUDE "iniprint.h"
    … …  
    158244    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1
    159245    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
     246    REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
    160247
    161248! In/Output arguments
    162249!****************************************************************************************
    163     REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
    164250    REAL, DIMENSION(klon), INTENT(INOUT) :: snow
    165     REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
    166251   
    167252! Output arguments
    … …  
    177262!****************************************************************************************
    178263    INTEGER               :: i,ki
     264    ! surface fluxes
    179265    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
    180     REAL, DIMENSION(klon) :: diff_sst, lmt_bils
     266    ! flux correction
     267    REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils
     268    ! surface wind stress
    181269    REAL, DIMENSION(klon) :: u0, v0
    182270    REAL, DIMENSION(klon) :: u1_lay, v1_lay
     271    ! ice creation
     272    REAL                  :: e_freeze, h_new, dfsic
    183273
    184274!****************************************************************************************
    … …  
    189279    beta(:)     = 1.
    190280    dif_grnd(:) = 0.
    191     agesno(:)   = 0.
    192281   
    193282! Suppose zero surface speed
    … …  
    215304    DO i=1,knon
    216305        ki=knindex(i)
    217         slab_bils(ki)=(fluxlat(i)+fluxsens(i)+radsol(i))*pctsrf(ki,is_oce)/(1.-zmasq(ki))
     306        slab_bils(ki)=(1.-fsic(ki))*(fluxlat(i)+fluxsens(i)+radsol(i) &
     307                      -precip_snow(i)*ice_lat*(1.+snow_wfact*fsic(ki)))
    218308        bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
    219309! Also taux, tauy, saved vars...
    … …  
    225315!****************************************************************************************
    226316    lmt_bils(:)=0.
    227     CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst) ! global pour un processus
    228     ! lmt_bils and diff_sst saved by limit_slab
    229     qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.
     317    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) ! global pour un processus
     318    ! lmt_bils and diff_sst,siv saved by limit_slab
     319    qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400.
    230320    ! qflux = total QFlux correction (in W/m2)
    231321
    … …  
    248338                tsurf_new(i)=tslab(ki,1)
    249339            END DO
    250         CASE('sicOBS') ! check for sea ice or tsurf below freezing
     340        CASE('sicOBS') ! check for sea ice or tslab below freezing
    251341            DO i=1,knon
    252342                ki=knindex(i)
    253                 IF ((tslab(ki,1).LT.t_freeze).OR.(pctsrf(ki,is_sic).GT.epsfra)) THEN
    254                     tsurf_new(i)=t_freeze
     343                IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
    255344                    tslab(ki,1)=t_freeze
    256                 ELSE
    257                     tsurf_new(i)=tslab(ki,1)
    258345                END IF
     346                tsurf_new(i)=tslab(ki,1)
    259347            END DO
    260348        CASE('sicINT')
    261349            DO i=1,knon
    262350                ki=knindex(i)
    263                 IF (pctsrf(ki,is_sic).LT.epsfra) THEN ! Free of ice
    264                     IF (tslab(ki,1).GT.t_freeze) THEN
     351                IF (fsic(ki).LT.epsfra) THEN ! Free of ice
     352                    IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
     353                        ! quantity of new ice formed
     354                        e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
     355                        ! new ice
     356                        tice(ki)=t_freeze
     357                        fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
     358                        IF (fsic(ki).GT.ice_frac_min) THEN
     359                            seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
     360                            tslab(ki,1)=t_freeze
     361                        ELSE
     362                            fsic(ki)=0.
     363                        END IF
     364                        tsurf_new(i)=t_freeze
     365                    ELSE
    265366                        tsurf_new(i)=tslab(ki,1)
    266                     ELSE
    267                         tsurf_new(i)=t_freeze
    268                         ! Call new ice routine
    269                         tslab(ki,1)=t_freeze
    270                     END IF
    271                 ELSE ! ice present, tslab update completed in slab_ice
     367                    END IF
     368                ELSE ! ice present
    272369                    tsurf_new(i)=t_freeze
    273                 END IF !ice free
     370                    IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
     371                        ! quantity of new ice formed over open ocean
     372                        e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
     373                                 /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
     374                        ! new ice height and fraction
     375                        h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
     376                        dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
     377                        h_new=MIN(e_freeze/dfsic,h_ice_max)
     378                        ! update tslab to freezing over open ocean only
     379                        tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki))
     380                        ! update sea ice
     381                        seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) &
     382                                   /(dfsic+fsic(ki))
     383                        fsic(ki)=fsic(ki)+dfsic
     384                        ! update snow?
     385                    END IF !freezing
     386                END IF ! sea ice present
    274387            END DO
    275388        END SELECT
    … …  
    279392!
    280393!****************************************************************************************
    281 !
    282 !  SUBROUTINE ocean_slab_ice(   &
    283 !       itime, dtime, jour, knon, knindex, &
    284 !       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
    285 !       AcoefH, AcoefQ, BcoefH, BcoefQ, &
    286 !       AcoefU, AcoefV, BcoefU, BcoefV, &
    287 !       ps, u1, v1, &
    288 !       radsol, snow, qsurf, qsol, agesno, tsoil, &
    289 !       alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
    290 !       tsurf_new, dflux_s, dflux_l)
    291 !
     394
     395  SUBROUTINE ocean_slab_ice(   &
     396       itime, dtime, jour, knon, knindex, &
     397       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
     398       AcoefH, AcoefQ, BcoefH, BcoefQ, &
     399       AcoefU, AcoefV, BcoefU, BcoefV, &
     400       ps, u1, v1, &
     401       radsol, snow, qsurf, qsol, agesno, &
     402       alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
     403       tsurf_new, dflux_s, dflux_l, swnet)
     404
     405   USE calcul_fluxs_mod
     406
     407   INCLUDE "YOMCST.h"
     408
     409! Input arguments
     410!****************************************************************************************
     411    INTEGER, INTENT(IN)                  :: itime, jour, knon
     412    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
     413    REAL, INTENT(IN)                     :: dtime
     414    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
     415    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
     416    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragm
     417    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
     418    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
     419    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
     420    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
     421    REAL, DIMENSION(klon), INTENT(IN)    :: ps
     422    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1
     423    REAL, DIMENSION(klon), INTENT(IN)    :: swnet
     424
     425! In/Output arguments
     426!****************************************************************************************
     427    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
     428    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
     429    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
     430
     431! Output arguments
     432!****************************************************************************************
     433    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
     434    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1_new  ! new albedo in visible SW interval
     435    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2_new  ! new albedo in near IR interval
     436    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
     437    REAL, DIMENSION(klon), INTENT(OUT)            :: flux_u1, flux_v1
     438    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
     439    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l
     440
     441! Local variables
     442!****************************************************************************************
     443    INTEGER               :: i,ki
     444    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
     445    REAL, DIMENSION(klon) :: u0, v0
     446    REAL, DIMENSION(klon) :: u1_lay, v1_lay
     447    ! intermediate heat fluxes:
     448    REAL                  :: f_cond, f_swpen
     449    ! for snow/ice albedo:
     450    REAL                  :: alb_snow, alb_ice, alb_pond
     451    REAL                  :: frac_snow, frac_ice, frac_pond
     452    ! for ice melt / freeze
     453    REAL                  :: e_melt, snow_evap, h_test
     454    ! dhsic, dfsic change in ice mass, fraction.
     455    REAL                  :: dhsic, dfsic, frac_mf
     456
    292457!****************************************************************************************
    293458! 1) Flux calculation
    294459!****************************************************************************************
    295 ! set beta, cal etc. depends snow / ice surf ?
     460! Suppose zero surface speed
     461    u0(:)=0.0
     462    v0(:)=0.0
     463    u1_lay(:) = u1(:) - u0(:)
     464    v1_lay(:) = v1(:) - v0(:)
     465
     466! set beta, cal, compute conduction fluxes inside ice/snow
     467    slab_bilg(:)=0.
     468    dif_grnd(:)=0.
     469    beta(:) = 1.
     470    DO i=1,knon
     471    ki=knindex(i)
     472        IF (snow(i).GT.snow_min) THEN
     473            ! snow-layer heat capacity
     474            cal(i)=2.*RCPD/(snow(i)*ice_cap)
     475            ! snow conductive flux
     476            f_cond=sno_cond*(tice(ki)-tsurf_in(i))/snow(i)
     477            ! all shortwave flux absorbed
     478            f_swpen=0.
     479            ! bottom flux (ice conduction)
     480            slab_bilg(ki)=ice_cond*(tice(ki)-t_freeze)/seaice(ki)
     481            ! update ice temperature
     482            tice(ki)=tice(ki)-2./ice_cap/(snow(i)+seaice(ki)) &
     483                     *(slab_bilg(ki)+f_cond)*dtime
     484       ELSE ! bare ice
     485            ! ice-layer heat capacity
     486            cal(i)=2.*RCPD/(seaice(ki)*ice_cap)
     487            ! conductive flux
     488            f_cond=ice_cond*(t_freeze-tice(ki))/seaice(ki)
     489            ! penetrative shortwave flux...
     490            f_swpen=swnet(i)*pen_frac*exp(-pen_ext*seaice(ki)/ice_den)
     491            slab_bilg(ki)=f_swpen-f_cond
     492        END IF
     493        radsol(i)=radsol(i)+f_cond-f_swpen
     494    END DO
     495    ! weight fluxes to ocean by sea ice fraction
     496    slab_bilg(:)=slab_bilg(:)*fsic(:)
     497
    296498! calcul_fluxs (sens, lat etc)
     499    CALL calcul_fluxs(knon, is_sic, dtime, &
     500        tsurf_in, p1lay, cal, beta, cdragh, ps, &
     501        precip_rain, precip_snow, snow, qsurf,  &
     502        radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
     503        AcoefH, AcoefQ, BcoefH, BcoefQ, &
     504        tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
     505    DO i=1,knon
     506        IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i)
     507    END DO
     508
    297509! calcul_flux_wind
    298 
    299 !****************************************************************************************
    300 ! 2) Update surface
    301 !****************************************************************************************
    302 ! neige, fonte
    303 ! flux glace-ocean
    304 ! update temperature
    305 ! neige precip, evap
    306 ! Melt snow & ice from above
     510    CALL calcul_flux_wind(knon, dtime, &
     511         u0, v0, u1, v1, cdragm, &
     512         AcoefU, AcoefV, BcoefU, BcoefV, &
     513         p1lay, temp_air, &
     514         flux_u1, flux_v1)
     515
     516!****************************************************************************************
     517! 2) Update snow and ice surface
     518!****************************************************************************************
     519! snow precip
     520    DO i=1,knon
     521        ki=knindex(i)
     522        IF (precip_snow(i) > 0.) THEN
     523            snow(i) = snow(i)+precip_snow(i)*dtime*(1.-snow_wfact*(1.-fsic(ki)))
     524        END IF
     525! snow and ice sublimation
     526        IF (evap(i) > 0.) THEN
     527           snow_evap = MIN (snow(i) / dtime, evap(i))
     528           snow(i) = snow(i) - snow_evap * dtime
     529           snow(i) = MAX(0.0, snow(i))
     530           seaice(ki) = MAX(0.0,seaice(ki)-(evap(i)-snow_evap)*dtime)
     531        ENDIF
     532! Melt / Freeze from above if Tsurf>0
     533        IF (tsurf_new(i).GT.t_melt) THEN
     534            ! energy available for melting snow (in kg/m2 of snow)
     535            e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. &
     536               /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
     537            ! remove snow
     538            IF (snow(i).GT.e_melt) THEN
     539                snow(i)=snow(i)-e_melt
     540                tsurf_new(i)=t_melt
     541            ELSE ! all snow is melted
     542                ! add remaining heat flux to ice
     543                e_melt=e_melt-snow(i)
     544                tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki))
     545                tsurf_new(i)=tice(ki)
     546            END IF
     547        END IF
     548! melt ice from above if Tice>0
     549        IF (tice(ki).GT.t_melt) THEN
     550            ! quantity of ice melted (kg/m2)
     551            e_melt=MAX(seaice(ki)*(tice(ki)-t_melt)*ice_cap/2. &
     552             /(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0)
     553            ! melt from above, height only
     554            dhsic=MIN(seaice(ki)-h_ice_min,e_melt)
     555            e_melt=e_melt-dhsic
     556            IF (e_melt.GT.0) THEN
     557            ! lateral melt if ice too thin
     558            dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/h_ice_min*fsic(ki))
     559            ! if all melted add remaining heat to ocean
     560            e_melt=MAX(0.,e_melt*fsic(ki)-dfsic*h_ice_min)
     561            slab_bilg(ki)=slab_bilg(ki)+ e_melt*ice_lat/dtime
     562            ! update height and fraction
     563            fsic(ki)=fsic(ki)-dfsic
     564            END IF
     565            seaice(ki)=seaice(ki)-dhsic
     566            ! surface temperature at melting point
     567            tice(ki)=t_melt
     568            tsurf_new(i)=t_melt
     569        END IF
     570        ! convert snow to ice if below floating line
     571        h_test=(seaice(ki)+snow(i))*ice_den-seaice(ki)*sea_den
     572        IF (h_test.GT.0.) THEN !snow under water
     573            ! extra snow converted to ice (with added frozen sea water)
     574            dhsic=h_test/(sea_den-ice_den+sno_den)
     575            seaice(ki)=seaice(ki)+dhsic
     576            snow(i)=snow(i)-dhsic*sno_den/ice_den
     577            ! available energy (freeze sea water + bring to tice)
     578            e_melt=dhsic*(1.-sno_den/ice_den)*(ice_lat+ &
     579                   ice_cap/2.*(t_freeze-tice(ki)))
     580            ! update ice temperature
     581            tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(i)+seaice(ki))
     582        END IF
     583    END DO
     584
    307585! New albedo
    308 
    309 !****************************************************************************************
    310 ! 3) Recalculate new ocean temperature
     586    DO i=1,knon
     587        ki=knindex(i)
     588       ! snow albedo: update snow age
     589        IF (snow(i).GT.0.0001) THEN
     590             agesno(i)=(agesno(i) + (1.-agesno(i)/50.)*dtime/86400.)&
     591                         * EXP(-1.*MAX(0.0,precip_snow(i))*dtime/5.)
     592        ELSE
     593            agesno(i)=0.0
     594        END IF
     595        ! snow albedo
     596        alb_snow=alb_sno_min+alb_sno_del*EXP(-agesno(i)/50.)
     597        ! ice albedo (varies with ice tkickness and temp)
     598        alb_ice=MAX(0.0,0.13*LOG(100.*seaice(ki)/ice_den)+0.1)
     599        IF (tice(ki).GT.t_freeze-0.01) THEN
     600            alb_ice=MIN(alb_ice,alb_ice_wet)
     601        ELSE
     602            alb_ice=MIN(alb_ice,alb_ice_dry)
     603        END IF
     604        ! pond albedo
     605        alb_pond=0.36-0.1*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
     606        ! pond fraction
     607        frac_pond=0.2*(2.0+MIN(0.0,MAX(tice(ki)-t_melt,-2.0)))
     608        ! snow fraction
     609        frac_snow=MAX(0.0,MIN(1.0-frac_pond,snow(i)/snow_min))
     610        ! ice fraction
     611        frac_ice=MAX(0.0,1.-frac_pond-frac_snow)
     612        ! total albedo
     613        alb1_new(i)=alb_snow*frac_snow+alb_ice*frac_ice+alb_pond*frac_pond
     614    END DO
     615    alb2_new(:) = alb1_new(:)
     616
     617!****************************************************************************************
     618! 3) Recalculate new ocean temperature (add fluxes below ice)
    311619!    Melt / freeze from below
    312620!***********************************************o*****************************************
    313 
    314 
    315 !  END SUBROUTINE ocean_slab_ice
     621    !cumul fluxes
     622    bilg_cum(:)=bilg_cum(:)+slab_bilg(:)
     623    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
     624        ! Add cumulated surface fluxes
     625        tslab(:,1)=tslab(:,1)+bilg_cum(:)*cyang*dtime
     626        DO i=1,knon
     627            ki=knindex(i)
     628            ! split lateral/top melt-freeze
     629            frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thin)/(h_ice_thick-h_ice_thin)))
     630            IF (tslab(ki,1).LE.t_freeze) THEN
     631               ! ****** Form new ice from below *******
     632               ! quantity of new ice
     633                e_melt=(t_freeze-tslab(ki,1))/cyang &
     634                       /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
     635               ! first increase height to h_thin
     636               dhsic=MAX(0.,MIN(h_ice_thin-seaice(ki),e_melt/fsic(ki)))
     637               seaice(ki)=dhsic+seaice(ki)
     638               e_melt=e_melt-fsic(ki)*dhsic
     639               IF (e_melt.GT.0.) THEN
     640               ! frac_mf fraction used for lateral increase
     641               dfsic=MIN(ice_frac_max-fsic(ki),e_melt*frac_mf/seaice(ki))
     642               fsic(ki)=fsic(ki)+dfsic
     643               e_melt=e_melt-dfsic*seaice(ki)
     644               ! rest used to increase height
     645               seaice(ki)=MIN(h_ice_max,seaice(ki)+e_melt/fsic(ki))
     646               END IF
     647               tslab(ki,1)=t_freeze
     648           ELSE ! slab temperature above freezing
     649               ! ****** melt ice from below *******
     650               ! quantity of melted ice
     651               e_melt=(tslab(ki,1)-t_freeze)/cyang &
     652                       /(ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
     653               ! first decrease height to h_thick
     654               dhsic=MAX(0.,MIN(seaice(ki)-h_ice_thick,e_melt/fsic(ki)))
     655               seaice(ki)=seaice(ki)-dhsic
     656               e_melt=e_melt-fsic(ki)*dhsic
     657               IF (e_melt.GT.0) THEN
     658               ! frac_mf fraction used for height decrease
     659               dhsic=MAX(0.,MIN(seaice(ki)-h_ice_min,e_melt*frac_mf/fsic(ki)))
     660               seaice(ki)=seaice(ki)-dhsic
     661               e_melt=e_melt-fsic(ki)*dhsic
     662               ! rest used to decrease fraction (up to 0!)
     663               dfsic=MIN(fsic(ki),e_melt/seaice(ki))
     664               ! keep remaining in ocean
     665               e_melt=e_melt-dfsic*seaice(ki)
     666               END IF
     667               tslab(ki,1)=t_freeze+e_melt*ice_lat*cyang
     668               fsic(ki)=fsic(ki)-dfsic
     669           END IF
     670        END DO
     671        bilg_cum(:)=0.
     672    END IF ! coupling time
     673   
     674    !tests fraction glace
     675    WHERE (fsic.LT.ice_frac_min)
     676        tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang
     677        tice=t_melt
     678        fsic=0.
     679        seaice=0.
     680    END WHERE
     681
     682  END SUBROUTINE ocean_slab_ice
    316683!
    317684!****************************************************************************************
    318685!
    319686  SUBROUTINE ocean_slab_final
    320   !, seaice_rst etc
    321 
    322     ! For ok_xxx vars (Ekman...)
    323     INCLUDE "clesphys.h"
    324687
    325688!****************************************************************************************
    326689! Deallocate module variables
    327 !
    328 !****************************************************************************************
    329     IF (ALLOCATED(pctsrf)) DEALLOCATE(pctsrf)
     690!****************************************************************************************
    330691    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
     692    IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
     693    IF (ALLOCATED(slab_bils)) DEALLOCATE(slab_bils)
     694    IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
     695    IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
     696    IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum)
     697    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
    331698
    332699  END SUBROUTINE ocean_slab_final

  • LMDZ5/branches/testing/libf/phylmd/ozonecm_m.F90

    r1999 r2220  
    9191    ozonecm = max(ozonecm, 1e-12)
    9292
    93     print*,'ozonecm Version2'
     93!    print*,'ozonecm Version2'
    9494
    9595  END function ozonecm

  • LMDZ5/branches/testing/libf/phylmd/pbl_surface_mod.F90

    r2187 r2220  
    3333  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE     :: fder   ! flux drift
    3434  !$OMP THREADPRIVATE(fder)
    35   REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: snow   ! snow at surface
     35  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: snow   ! snow at surface
    3636  !$OMP THREADPRIVATE(snow)
    3737  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: qsurf  ! humidity at surface
    … …  
    172172       debut,     lafin,                              &
    173173       rlon,      rlat,      rugoro,   rmu0,          &
    174        zsig,      sollwd_m,  pphi,     cldt,          &
     174       zsig,      lwdown_m,  pphi,     cldt,          &
    175175       rain_f,    snow_f,    solsw_m,  sollw_m,       &
    176176       t,         q,         u,        v,             &
    … …  
    182182       pplay,     paprs,     pctsrf,                  &
    183183       ts,        alb1, alb2,ustar, u10m, v10m,wstar, &
    184        lwdown_m,  cdragh,    cdragm,   zu1,    zv1,   &
     184       cdragh,    cdragm,   zu1,    zv1,              &
    185185       alb1_m,    alb2_m,    zxsens,   zxevap,        &
    186186       alb3_lic,  runoff,    snowhgt,   qsnow,     to_ice,    sissnow,  &
    … …  
    327327! Martin
    328328    REAL, DIMENSION(klon),        INTENT(IN)        :: zsig    ! slope
    329     REAL, DIMENSION(klon),        INTENT(IN)        :: sollwd_m ! net longwave radiation at mean s   
     329    REAL, DIMENSION(klon),        INTENT(IN)        :: lwdown_m ! downward longwave radiation at mean s   
    330330    REAL, DIMENSION(klon),        INTENT(IN)        :: cldt    ! total cloud fraction
    331331    REAL, DIMENSION(klon,klev),   INTENT(IN)        :: pphi    ! geopotential (m2/s2)
    … …  
    367367! Output variables
    368368!****************************************************************************************
    369     REAL, DIMENSION(klon),        INTENT(OUT)       :: lwdown_m   ! Downcoming longwave radiation
    370369    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragh     ! drag coefficient for T and Q
    371370    REAL, DIMENSION(klon),        INTENT(OUT)       :: cdragm     ! drag coefficient for wind
    … …  
    780779    ! Martin
    781780    REAL, DIMENSION(klon, nbsrf)       :: sollwd ! net longwave radiation at surface
    782     REAL, DIMENSION(klon)              :: ysollwd
    783781    REAL, DIMENSION(klon)              :: ytoice
    784782    REAL, DIMENSION(klon)              :: ysnowhgt, yqsnow, ysissnow, yrunoff
    … …  
    855853! 2a) Initialization of all argument variables with INTENT(OUT)
    856854!****************************************************************************************
    857  lwdown_m(:)=0.
    858855 cdragh(:)=0. ; cdragm(:)=0.
    859856 zu1(:)=0. ; zv1(:)=0.
    … …  
    938935    ! Martin
    939936    ysnowhgt = 0.0; yqsnow = 0.0     ; yrunoff = 0.0   ; ytoice =0.0
    940     yalb3_new = 0.0  ; ysissnow = 0.0  ; ysollwd = 0.0
     937    yalb3_new = 0.0  ; ysissnow = 0.0
    941938    ypphi = 0.0   ; ycldt = 0.0      ; yrmu0 = 0.0
    942939    ! Martin
    … …  
    11091106       DO i = 1, klon
    11101107          sollw(i,nsrf) = sollw_m(i) + 4.0*RSIGMA*ztsol(i)**3 * (ztsol(i)-ts(i,nsrf))
    1111           ! Martin
    1112           sollwd(i,nsrf)= sollwd_m(i)
    1113           ! Martin
     1108
     1109!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1110!         ! Martin
     1111! Apparently introduced for sisvat but not used
     1112!         sollwd(i,nsrf)= sollwd_m(i)
     1113!         ! Martin
     1114!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1115
    11141116          solsw(i,nsrf) = solsw_m(i) * (1.-alb(i,nsrf)) / (1.-alb_m(i))
    11151117       ENDDO
    11161118    ENDDO
    11171119
    1118 
    1119 ! Downwelling longwave radiation at mean surface
    1120     lwdown_m(:) = 0.0
    1121     DO i = 1, klon
    1122        lwdown_m(i) = sollw_m(i) + RSIGMA*ztsol(i)**4
    1123     ENDDO
    11241120
    11251121!****************************************************************************************
    … …  
    11801176          yagesno(j) = agesno(i,nsrf)
    11811177          yfder(j)   = fder(i)
     1178          ylwdown(j) = lwdown_m(i)
    11821179          ysolsw(j)  = solsw(i,nsrf)
    11831180          ysollw(j)  = sollw(i,nsrf)
    … …  
    17031700     
    17041701       CASE(is_ter)
    1705           ! ylwdown : to be removed, calculation is now done at land surface in surf_land
    1706           ylwdown(:)=0.0
    1707           DO i=1,knon
    1708              ylwdown(i)=lwdown_m(ni(i))
    1709           END DO
     1702!          print*,"DEBUGTS",yts(knon/2),ylwdown(knon/2)
    17101703          CALL surf_land(itap, dtime, date0, jour, knon, ni,&
    17111704               rlon, rlat, &
    … …  
    17461739          CALL surf_landice(itap, dtime, knon, ni, &
    17471740               rlon, rlat, debut, lafin, &
    1748                yrmu0, ysollwd, yalb, ypphi(:,1), &
     1741               yrmu0, ylwdown, yalb, ypphi(:,1), &
    17491742               ysolsw, ysollw, yts, ypplay(:,1), &
    17501743               ycdragh, ycdragm, yrain_f, ysnow_f, yt(:,1), yq(:,1),&

  • LMDZ5/branches/testing/libf/phylmd/phyaqua_mod.F90

    r2160 r2220  
    238238    rugos = rugos_omp
    239239    WRITE (*, *) 'iniaqua: rugos=', rugos
    240     zmasq(:) = pctsrf(:, is_oce)
     240    zmasq(:) = pctsrf(:, is_ter)
    241241
    242242    ! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
    … …  
    538538      dims(2) = ntim
    539539
    540       ! cc      ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
     540#ifdef NC_DOUBLE
     541      ierr = nf_def_var(nid, 'TEMPS', nf_double, 1, ntim, id_tim)
     542#else
    541543      ierr = nf_def_var(nid, 'TEMPS', nf_float, 1, ntim, id_tim)
     544#endif
    542545      ierr = nf_put_att_text(nid, id_tim, 'title', 17, 'Jour dans l annee')
    543       ! cc      ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
     546
     547#ifdef NC_DOUBLE
     548      ierr = nf_def_var(nid, 'NAT', nf_double, 2, dims, id_nat)
     549#else
    544550      ierr = nf_def_var(nid, 'NAT', nf_float, 2, dims, id_nat)
     551#endif
    545552      ierr = nf_put_att_text(nid, id_nat, 'title', 23, &
    546553        'Nature du sol (0,1,2,3)')
    547       ! cc      ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
     554
     555#ifdef NC_DOUBLE
     556      ierr = nf_def_var(nid, 'SST', nf_double, 2, dims, id_sst)
     557#else
    548558      ierr = nf_def_var(nid, 'SST', nf_float, 2, dims, id_sst)
     559#endif
    549560      ierr = nf_put_att_text(nid, id_sst, 'title', 35, &
    550561        'Temperature superficielle de la mer')
    551       ! cc      ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
     562
     563#ifdef NC_DOUBLE
     564      ierr = nf_def_var(nid, 'BILS', nf_double, 2, dims, id_bils)
     565#else
    552566      ierr = nf_def_var(nid, 'BILS', nf_float, 2, dims, id_bils)
     567#endif
    553568      ierr = nf_put_att_text(nid, id_bils, 'title', 32, &
    554569        'Reference flux de chaleur au sol')
    555       ! cc      ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
     570
     571#ifdef NC_DOUBLE
     572      ierr = nf_def_var(nid, 'ALB', nf_double, 2, dims, id_alb)
     573#else
    556574      ierr = nf_def_var(nid, 'ALB', nf_float, 2, dims, id_alb)
     575#endif
    557576      ierr = nf_put_att_text(nid, id_alb, 'title', 19, 'Albedo a la surface')
    558       ! cc      ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
     577
     578#ifdef NC_DOUBLE
     579      ierr = nf_def_var(nid, 'RUG', nf_double, 2, dims, id_rug)
     580#else
    559581      ierr = nf_def_var(nid, 'RUG', nf_float, 2, dims, id_rug)
     582#endif
    560583      ierr = nf_put_att_text(nid, id_rug, 'title', 8, 'Rugosite')
    561584
     585#ifdef NC_DOUBLE
     586      ierr = nf_def_var(nid, 'FTER', nf_double, 2, dims, id_fter)
     587#else
    562588      ierr = nf_def_var(nid, 'FTER', nf_float, 2, dims, id_fter)
    563       ierr = nf_put_att_text(nid, id_fter, 'title', 8, 'Frac. Terre')
     589#endif
     590      ierr = nf_put_att_text(nid, id_fter, 'title',10,'Frac. Land')
     591#ifdef NC_DOUBLE
     592      ierr = nf_def_var(nid, 'FOCE', nf_double, 2, dims, id_foce)
     593#else
    564594      ierr = nf_def_var(nid, 'FOCE', nf_float, 2, dims, id_foce)
    565       ierr = nf_put_att_text(nid, id_foce, 'title', 8, 'Frac. Terre')
     595#endif
     596      ierr = nf_put_att_text(nid, id_foce, 'title',11,'Frac. Ocean')
     597#ifdef NC_DOUBLE
     598      ierr = nf_def_var(nid, 'FSIC', nf_double, 2, dims, id_fsic)
     599#else
    566600      ierr = nf_def_var(nid, 'FSIC', nf_float, 2, dims, id_fsic)
    567       ierr = nf_put_att_text(nid, id_fsic, 'title', 8, 'Frac. Terre')
     601#endif
     602      ierr = nf_put_att_text(nid, id_fsic, 'title',13,'Frac. Sea Ice')
     603#ifdef NC_DOUBLE
     604      ierr = nf_def_var(nid, 'FLIC', nf_double, 2, dims, id_flic)
     605#else
    568606      ierr = nf_def_var(nid, 'FLIC', nf_float, 2, dims, id_flic)
    569       ierr = nf_put_att_text(nid, id_flic, 'title', 8, 'Frac. Terre')
     607#endif
     608      ierr = nf_put_att_text(nid, id_flic, 'title',14,'Frac. Land Ice')
    570609
    571610      ierr = nf_enddef(nid)

  • LMDZ5/branches/testing/libf/phylmd/phyetat0.F90

    r2187 r2220  
    88  USE fonte_neige_mod,  ONLY : fonte_neige_init
    99  USE pbl_surface_mod,  ONLY : pbl_surface_init
    10   USE surface_data,     ONLY : type_ocean
     10  USE surface_data,     ONLY : type_ocean, version_ocean
    1111  USE phys_state_var_mod, ONLY : ancien_ok, clwcon, detr_therm, dtime, &
    1212       du_gwd_rando, dv_gwd_rando, entr_therm, f0, falb1, falb2, fm_therm, &
    1313       ftsol, pbl_tke, pctsrf, q_ancien, radpas, radsol, rain_fall, ratqs, &
    14        rlat, rlon, rnebcon, rugoro, sig1, snow_fall, solaire_etat0, sollw, &
     14       rlat, rlon, rnebcon, rugoro, sig1, snow_fall, solaire_etat0, sollw, sollwdown, &
    1515       solsw, t_ancien, u_ancien, v_ancien, w01, wake_cstar, wake_deltaq, &
    1616       wake_deltat, wake_delta_pbl_TKE, delta_tsurf, wake_fip, wake_pe, &
    … …  
    2323  USE carbon_cycle_mod, ONLY : carbon_cycle_tr, carbon_cycle_cpl, co2_send
    2424  USE indice_sol_mod, only: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
    25   USE ocean_slab_mod, ONLY: tslab, ocean_slab_init
     25  USE ocean_slab_mod, ONLY: tslab, seaice, tice, ocean_slab_init
    2626
    2727  IMPLICIT none
    … …  
    3737  include "thermcell.h"
    3838  include "compbl.h"
     39  include "YOMCST.h"
    3940  !======================================================================
    4041  CHARACTER*(*) fichnom
    … …  
    5354  REAL fractint(klon)
    5455  REAL trs(klon, nbtr)
     56  REAL zts(klon)
    5557
    5658  CHARACTER*6 ocean_in
    … …  
    513515  PRINT*, 'Rayonnement IF au sol sollw:', xmin, xmax
    514516
     517  CALL get_field("sollwdown", sollwdown, found)
     518  IF (.NOT. found) THEN
     519     PRINT*, 'phyetat0: Le champ <sollwdown> est absent'
     520     PRINT*, 'mis a zero'
     521     sollwdown = 0.
     522     zts=0.
     523     do nsrf=1,nbsrf
     524        zts(:)=zts(:)+ftsol(:,nsrf)*pctsrf(:,nsrf)
     525     enddo
     526     sollwdown(:)=sollw(:)+RSIGMA*zts(:)**4
     527  ENDIF
     528!  print*,'TS SOLL',zts(klon/2),sollw(klon/2),sollwdown(klon/2)
     529  xmin = 1.0E+20
     530  xmax = -1.0E+20
     531  DO i = 1, klon
     532     xmin = MIN(sollwdown(i), xmin)
     533     xmax = MAX(sollwdown(i), xmax)
     534  ENDDO
     535  PRINT*, 'Rayonnement IF au sol sollwdown:', xmin, xmax
     536
     537
    515538  ! Lecture derive des flux:
    516539
    … …  
    11371160  ! Initialize Slab variables
    11381161  IF ( type_ocean == 'slab' ) THEN
    1139       ALLOCATE(tslab(klon,nslay), stat=ierr)
    1140       IF (ierr /= 0) CALL abort_gcm &
    1141          ('phyetat0', 'pb allocation tslab', 1)
     1162      print*, "calling slab_init"
     1163      CALL ocean_slab_init(dtime, pctsrf)
     1164      ! tslab
    11421165      CALL get_field("tslab", tslab, found)
    11431166      IF (.NOT. found) THEN
    … …  
    11451168          PRINT*, "Initialisation a tsol_oce"
    11461169          DO i=1,nslay
    1147               tslab(:,i)=ftsol(:,is_oce)
     1170              tslab(:,i)=MAX(ftsol(:,is_oce),271.35)
    11481171          END DO
    11491172      END IF
    1150       print*, "calling slab_init"
    1151       CALL ocean_slab_init(dtime, pctsrf)
     1173      ! Sea ice variables
     1174      IF (version_ocean == 'sicINT') THEN
     1175          CALL get_field("slab_tice", tice, found)
     1176          IF (.NOT. found) THEN
     1177              PRINT*, "phyetat0: Le champ <tice> est absent"
     1178              PRINT*, "Initialisation a tsol_sic"
     1179                  tice(:)=ftsol(:,is_sic)
     1180          END IF
     1181          CALL get_field("seaice", seaice, found)
     1182          IF (.NOT. found) THEN
     1183              PRINT*, "phyetat0: Le champ <seaice> est absent"
     1184              PRINT*, "Initialisation a 0/1m suivant fraction glace"
     1185              seaice(:)=0.
     1186              WHERE (pctsrf(:,is_sic).GT.EPSFRA)
     1187                  seaice=917.
     1188              END WHERE
     1189          END IF
     1190      END IF !sea ice INT
    11521191  END IF ! Slab       
    11531192

  • LMDZ5/branches/testing/libf/phylmd/phyredem.F90

    r2073 r2220  
    1616  USE indice_sol_mod
    1717  USE surface_data
    18   USE ocean_slab_mod, ONLY : tslab
     18  USE ocean_slab_mod, ONLY : tslab, seaice, tice, fsic
    1919
    2020  IMPLICIT none
    … …  
    107107  ! BP ajout des fraction de chaque sous-surface
    108108
     109  ! Get last fractions from slab ocean
     110  IF (type_ocean == 'slab' .AND. version_ocean == "sicINT") THEN
     111      WHERE (1.-zmasq(:).GT.EPSFRA)
     112          pctsrf(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:))
     113          pctsrf(:,is_sic)=fsic(:)*(1.-zmasq(:))
     114      END WHERE
     115  END IF
     116
    109117  ! 1. fraction de terre
    110118
    … …  
    209217
    210218  CALL put_field("sollw", "Rayonnement IF a la surface", sollw)
     219
     220  CALL put_field("sollwdown", "Rayonnement down IF a la surface", sollw)
    211221
    212222  CALL put_field("fder", "Derive de flux", fder)
    … …  
    350360  IF (type_ocean == 'slab') THEN
    351361      CALL put_field("tslab", "Slab ocean temperature", tslab)
     362      IF (version_ocean == 'sicINT') THEN
     363          CALL put_field("seaice", "Slab seaice (kg/m2)", seaice)
     364          CALL put_field("slab_tice", "Slab sea ice temperature", tice)
     365      END IF
    352366  END IF
    353367

  • LMDZ5/branches/testing/libf/phylmd/phys_local_var_mod.F90

    r2187 r2220  
    3535      REAL, SAVE, ALLOCATABLE :: d_t_lsc(:,:),d_q_lsc(:,:),d_ql_lsc(:,:),d_qi_lsc(:,:)
    3636      !$OMP THREADPRIVATE(d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc)
     37      REAL, SAVE, ALLOCATABLE :: d_t_lwr(:,:),d_t_lw0(:,:),d_t_swr(:,:),d_t_sw0(:,:)
     38      !$OMP THREADPRIVATE(d_t_lwr,d_t_lw0,d_t_swr,d_t_sw0)
    3739      REAL, SAVE, ALLOCATABLE :: d_t_ajsb(:,:), d_q_ajsb(:,:)
    3840      !$OMP THREADPRIVATE(d_t_ajsb, d_q_ajsb)
    … …  
    393395      allocate(d_t_wake(klon,klev),d_q_wake(klon,klev))
    394396      allocate(d_t_lsc(klon,klev),d_q_lsc(klon,klev))
     397      allocate(d_t_lwr(klon,klev),d_t_lw0(klon,klev))
     398      allocate(d_t_swr(klon,klev),d_t_sw0(klon,klev))
    395399      allocate(d_ql_lsc(klon,klev),d_qi_lsc(klon,klev))
    396400      allocate(d_t_ajsb(klon,klev),d_q_ajsb(klon,klev))
    … …  
    599603      deallocate(d_t_wake,d_q_wake)
    600604      deallocate(d_t_lsc,d_q_lsc)
     605      deallocate(d_t_lwr,d_t_lw0)
     606      deallocate(d_t_swr,d_t_sw0)
    601607      deallocate(d_ql_lsc,d_qi_lsc)
    602608      deallocate(d_t_ajsb,d_q_ajsb)

  • LMDZ5/branches/testing/libf/phylmd/phys_output_ctrlout_mod.F90

    r2187 r2220  
    465465  TYPE(ctrl_out), SAVE :: o_slab_bils = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
    466466    'slab_bils', 'flux atmos - slab ponderes foce', 'W/m2', (/ ('', i=1, 9) /))
     467  TYPE(ctrl_out), SAVE :: o_slab_bilg = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
     468    'slab_bilg', 'flux glace - slab ponderes fsic', 'W/m2', (/ ('', i=1, 9) /))
    467469  TYPE(ctrl_out), SAVE :: o_slab_qflux = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
    468470    'slab_qflux', 'Correction flux slab', 'W/m2', (/ ('', i=1, 9) /))
    469471  TYPE(ctrl_out), SAVE :: o_tslab = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
    470472    'tslab', 'Temperature ocean slab', 'K', (/ ('', i=1, 9) /))
     473  TYPE(ctrl_out), SAVE :: o_slab_tice = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
     474    'slab_tice', 'Temperature banquise slab', 'K', (/ ('', i=1, 9) /))
     475  TYPE(ctrl_out), SAVE :: o_slab_sic = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
     476    'seaice', 'Epaisseur banquise slab', 'kg/m2', (/ ('', i=1, 9) /))
    471477  TYPE(ctrl_out), SAVE :: o_ale_bl = ctrl_out((/ 1, 1, 1, 10, 10, 10, 11, 11, 11 /), &
    472478    'ale_bl', 'ALE BL', 'm2/s2', (/ ('', i=1, 9) /))
    … …  
    10081014      ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11 /),'ages_oce',"Snow age", "day", (/ ('', i=1, 9) /)), &
    10091015      ctrl_out((/ 3, 10, 10, 10, 10, 10, 11, 11, 11 /),'ages_sic',"Snow age", "day", (/ ('', i=1, 9) /)) /)
     1016
     1017  TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_snow_srf     = (/ &
     1018      ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11 /),'snow_ter', "Snow", "kg/m2", (/ ('', i=1, 9) /)), &
     1019      ctrl_out((/ 3, 10, 10, 10, 10, 10, 11, 11, 11 /),'snow_lic', "Snow", "kg/m2", (/ ('', i=1, 9) /)), &
     1020      ctrl_out((/ 10, 10, 10, 10, 10, 10, 11, 11, 11 /),'snow_oce',"Snow", "kg/m2", (/ ('', i=1, 9) /)), &
     1021      ctrl_out((/ 3, 10, 10, 10, 10, 10, 11, 11, 11 /),'snow_sic',"Snow", "kg/m2", (/ ('', i=1, 9) /)) /)
    10101022
    10111023  TYPE(ctrl_out), SAVE, DIMENSION(4) :: o_rugs_srf     = (/ &

  • LMDZ5/branches/testing/libf/phylmd/phys_output_write_mod.F90

    r2187 r2220  
    8080         o_alp_bl_conv, o_alp_bl_stat, &
    8181         o_slab_qflux, o_tslab, o_slab_bils, &
     82         o_slab_bilg, o_slab_sic, o_slab_tice, &
    8283         o_weakinv, o_dthmin, o_cldtau, &
    8384         o_cldemi, o_pr_con_l, o_pr_con_i, &
    … …  
    113114         o_rnebls, o_rhum, o_ozone, o_ozone_light, &
    114115         o_dtphy, o_dqphy, o_albe_srf, o_rugs_srf, &
    115          o_ages_srf, o_alb1, o_alb2, o_tke, &
     116         o_ages_srf, o_snow_srf, o_alb1, o_alb2, o_tke, &
    116117         o_tke_max, o_kz, o_kz_max, o_clwcon, &
    117118         o_dtdyn, o_dqdyn, o_dudyn, o_dvdyn, &
    … …  
    166167         wake_deltaq, ftd, fqd, ale_bl_trig, albsol1, &
    167168         rnebcon, wo, falb1, albsol2, coefh, clwcon0, &
    168          ratqs, entr_therm, zqasc, detr_therm, f0, heat, &
    169          heat0, cool, cool0, lwup, lwdn, lwup0, coefm, &
     169         ratqs, entr_therm, zqasc, detr_therm, f0, &
     170         lwup, lwdn, lwup0, coefm, &
    170171         swupp, lwupp, swup0p, lwup0p, swdnp, lwdnp, &
    171172         swdn0p, lwdn0p, tnondef, O3sumSTD, uvsumSTD, &
    … …  
    215216         d_u_ajs, d_v_ajs, &
    216217         d_u_con, d_v_con, d_q_con, d_q_ajs, d_t_lsc, &
     218         d_t_lwr,d_t_lw0,d_t_swr,d_t_sw0, &
    217219         d_t_eva, d_q_lsc, beta_prec, d_t_lscth, &
    218220         d_t_lscst, d_q_lscth, d_q_lscst, plul_th, &
    … …  
    226228         bils_ec,bils_ech, bils_tke, bils_kinetic, bils_latent, bils_enthalp, &
    227229         itau_con, nfiles, clef_files, nid_files, zvstr_gwd_rando
    228     USE ocean_slab_mod, only: tslab, slab_bils
     230    USE ocean_slab_mod, only: tslab, slab_bils, slab_bilg, tice, seaice
     231    USE pbl_surface_mod, only: snow
    229232    USE indice_sol_mod, only: nbsrf
    230233    USE infotrac, only: nqtot, nqo, type_trac
    231234    USE comgeomphy, only: airephy
    232     USE surface_data, only: type_ocean, ok_veget, ok_snow
     235    USE surface_data, only: type_ocean, version_ocean, ok_veget, ok_snow
    233236!    USE aero_mod, only: naero_spc
    234237    USE aero_mod, only: naero_tot, id_STRAT_phy
    … …  
    399402       CALL histwrite_phy(o_pluc, zx_tmp_fi2d)
    400403       CALL histwrite_phy(o_snow, snow_fall)
    401        CALL histwrite_phy(o_msnow, snow_o)
     404       CALL histwrite_phy(o_msnow, zxsnow)
    402405       CALL histwrite_phy(o_fsnow, zfra_o)
    403406       CALL histwrite_phy(o_evap, evap)
    … …  
    513516
    514517       IF (ok_snow) THEN
    515           CALL histwrite_phy(o_snowsrf, zxsnow)
     518          CALL histwrite_phy(o_snowsrf, snow_o)
    516519          CALL histwrite_phy(o_qsnow, qsnow)
    517520          CALL histwrite_phy(o_snowhgt,snowhgt)
    … …  
    754757              CALL histwrite_phy(o_tslab, tslab)
    755758          END IF
     759          IF (version_ocean=='sicINT') THEN
     760              CALL histwrite_phy(o_slab_bilg, slab_bilg)
     761              CALL histwrite_phy(o_slab_tice, tice)
     762              CALL histwrite_phy(o_slab_sic, seaice)
     763          END IF
    756764       ENDIF !type_ocean == force/slab
    757765       CALL histwrite_phy(o_weakinv, weak_inversion)
    … …  
    969977          IF (vars_defined) zx_tmp_fi2d(1 : klon) = agesno( 1 : klon, nsrf)
    970978          CALL histwrite_phy(o_ages_srf(nsrf), zx_tmp_fi2d)
     979          IF (vars_defined) zx_tmp_fi2d(1 : klon) = snow( 1 : klon, nsrf)
     980          CALL histwrite_phy(o_snow_srf(nsrf), zx_tmp_fi2d)
    971981       ENDDO !nsrf=1, nbsrf
    972982       CALL histwrite_phy(o_alb1, albsol1)
    … …  
    11271137       IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_q_ajsb(1:klon,1:klev)/pdtphys
    11281138       CALL histwrite_phy(o_dqajs, zx_tmp_fi3d)
    1129        IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=heat(1:klon,1:klev)/RDAY
     1139       IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_swr(1:klon,1:klev)/pdtphys
    11301140       CALL histwrite_phy(o_dtswr, zx_tmp_fi3d)
    1131        IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=heat0(1:klon,1:klev)/RDAY
     1141       IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_sw0(1:klon,1:klev)/pdtphys
    11321142       CALL histwrite_phy(o_dtsw0, zx_tmp_fi3d)
    1133        IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=-1.*cool(1:klon,1:klev)/RDAY
     1143       IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_lwr(1:klon,1:klev)/pdtphys
    11341144       CALL histwrite_phy(o_dtlwr, zx_tmp_fi3d)
    1135        IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=-1.*cool0(1:klon,1:klev)/RDAY
     1145       IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_lw0(1:klon,1:klev)/pdtphys
    11361146       CALL histwrite_phy(o_dtlw0, zx_tmp_fi3d)
    11371147       IF(vars_defined) zx_tmp_fi3d(1:klon,1:klev)=d_t_ec(1:klon,1:klev)/pdtphys
    … …  
    11931203       CALL histwrite_phy(o_tnt, zx_tmp_fi3d)
    11941204       IF(vars_defined) THEN
    1195           zx_tmp_fi3d(1:klon,1:klev)=heat(1:klon,1:klev)/RDAY - &
    1196                cool(1:klon,1:klev)/RDAY
     1205          zx_tmp_fi3d(1:klon,1:klev)=d_t_swr(1:klon,1:klev)/pdtphys + &
     1206               d_t_lwr(1:klon,1:klev)/pdtphys
    11971207       ENDIF
    11981208       CALL histwrite_phy(o_tntr, zx_tmp_fi3d)

  • LMDZ5/branches/testing/libf/phylmd/phys_state_var_mod.F90

    r2187 r2220  
    6060      REAL, ALLOCATABLE, SAVE :: clwcon(:,:),rnebcon(:,:)
    6161!$OMP THREADPRIVATE(clwcon,rnebcon)
     62      REAL, ALLOCATABLE, SAVE :: qtc_cv(:,:),sigt_cv(:,:)
     63!$OMP THREADPRIVATE(qtc_cv,sigt_cv)
    6264      REAL, ALLOCATABLE, SAVE :: ratqs(:,:)
    6365!$OMP THREADPRIVATE(ratqs)
    … …  
    416418!!! Rom P <<<
    417419      ALLOCATE(clwcon(klon,klev),rnebcon(klon,klev))
     420      ALLOCATE(qtc_cv(klon,klev),sigt_cv(klon,klev))
    418421      ALLOCATE(ratqs(klon,klev))
    419422      ALLOCATE(pbl_tke(klon,klev+1,nbsrf+1))
    … …  
    566569      deallocate(zthe, zpic, zval)
    567570      deallocate(rugoro, t_ancien, q_ancien, clwcon, rnebcon)
     571      deallocate(qtc_cv,sigt_cv)
    568572      deallocate(        u_ancien, v_ancien                 )
    569573      deallocate(        tr_ancien)                           !RomP

  • LMDZ5/branches/testing/libf/phylmd/physiq.F90

    r2187 r2220  
    636636  !$OMP THREADPRIVATE(fact_cldcon,facttemps)
    637637
    638   integer iflag_cldcon
    639   save iflag_cldcon
    640   !$OMP THREADPRIVATE(iflag_cldcon)
     638  integer iflag_cldth
     639  save iflag_cldth
     640  !$OMP THREADPRIVATE(iflag_cldth)
    641641  logical ptconv(klon,klev)
    642642  !IM cf. AM 081204 BEG
    … …  
    913913          solarlong0,seuil_inversion, &
    914914          fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
    915           iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
     915          iflag_cldth,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
    916916          ok_ade, ok_aie, ok_cdnc, aerosol_couple,  &
    917917          flag_aerosol, flag_aerosol_strat, new_aod, &
    … …  
    10141014     print*,'CYCLE_DIURNE', cycle_diurne
    10151015     !
    1016      IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
    1017         abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
     1016     IF (iflag_con.EQ.2.AND.iflag_cldth.GT.-1) THEN
     1017        abort_message = 'Tiedtke needs iflag_cldth=-2 or -1'
    10181018        CALL abort_gcm (modname,abort_message,1)
    10191019     ENDIF
    … …  
    11301130                ,alp_bl_prescr, ale_bl_prescr)
    11311131           ! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
    1132            !        print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
     1132           !        print*,'apres ini_wake iflag_cldth=', iflag_cldth
    11331133        endif
    11341134
    … …  
    13081308             annee_ref, &
    13091309             day_ref,  &
     1310             day_ini, &
     1311             start_time, &
    13101312             itau_phy, &
    13111313             io_lon, &
    … …  
    18111813!>nrlmd+jyg
    18121814          pplay,     paprs,     pctsrf,             &
    1813           ftsol,falb1,falb2,ustar,u10m,v10m,wstar, &
    1814           sollwdown, cdragh,    cdragm,  u1,    v1, &
    1815           albsol1,   albsol2,   sens,    evap,   &
     1815          ftsol,falb1,falb2,ustar,u10m,v10m,wstar,  &
     1816          cdragh,    cdragm,  u1,    v1,            &
     1817          albsol1,   albsol2,   sens,    evap,      &
    18161818          albsol3_lic,runoff,   snowhgt,   qsnow, to_ice, sissnow, &
    18171819          zxtsol,    zxfluxlat, zt2m,    qsat2m,  &
    … …  
    21682170             ftd,fqd,lalim_conv,wght_th, &
    21692171             ev, ep,epmlmMm,eplaMm, &
    2170              wdtrainA,wdtrainM,wght_cvfd)
     2172             wdtrainA,wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
     2173             tau_cld_cv,coefw_cld_cv)
    21712174        ! RomP <<<
    21722175
    … …  
    22182221     !   calcul des proprietes des nuages convectifs
    22192222     clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
     2223     IF (iflag_cld_cv <= 1) THEN
    22202224     call clouds_gno &
    22212225          (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
     2226     ELSE
     2227     call clouds_bigauss &
     2228          (klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
     2229     ENDIF
     2230
    22222231
    22232232     ! =================================================================== c
    … …  
    24522461  END IF
    24532462
    2454   !      print*,'apres callwake iflag_cldcon=', iflag_cldcon
     2463  !      print*,'apres callwake iflag_cldth=', iflag_cldth
    24552464  !
    24562465  !===================================================================
    … …  
    26222631              enddo
    26232632
    2624            ELSE IF (iflag_trig_bl.eq.2) then
     2633           ELSE IF (iflag_trig_bl.ge.2) then
    26252634
    26262635              do i=1,klon
    … …  
    27732782  ! water distribution
    27742783  CALL  calcratqs(klon,klev,prt_level,lunout,        &
    2775        iflag_ratqs,iflag_con,iflag_cldcon,pdtphys,  &
     2784       iflag_ratqs,iflag_con,iflag_cldth,pdtphys,  &
    27762785       ratqsbas,ratqshaut,tau_ratqs,fact_cldcon,   &
    27772786       ptconv,ptconvth,clwcon0th, rnebcon0th,     &
    … …  
    27952804       frac_impa, frac_nucl, beta_prec_fisrt, &
    27962805       prfl, psfl, rhcl,  &
    2797        zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon, &
     2806       zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldth, &
    27982807       iflag_ice_thermo)
    27992808  !
    … …  
    28512860  !
    28522861  !IM cf FH
    2853   !     IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
    2854   IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
     2862  !     IF (iflag_cldth.eq.-1) THEN ! seulement pour Tiedtke
     2863  IF (iflag_cldth.le.-1) THEN ! seulement pour Tiedtke
    28552864     snow_tiedtke=0.
    28562865     !     print*,'avant calcul de la pseudo precip '
    2857      !     print*,'iflag_cldcon',iflag_cldcon
    2858      if (iflag_cldcon.eq.-1) then
     2866     !     print*,'iflag_cldth',iflag_cldth
     2867     if (iflag_cldth.eq.-1) then
    28592868        rain_tiedtke=rain_con
    28602869     else
    … …  
    28892898     ENDDO
    28902899
    2891   ELSE IF (iflag_cldcon.ge.3) THEN
     2900  ELSE IF (iflag_cldth.ge.3) THEN
    28922901     !  On prend pour les nuages convectifs le max du calcul de la
    28932902     !  convection et du calcul du pas de temps precedent diminue d'un facteur
    … …  
    29322941             tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,  &
    29332942             tausum_aero, tau3d_aero)
     2943
     2944             CALL aeropt_lw_rrtm
    29342945#else
    29352946
    … …  
    29762987     !   On prend la somme des fractions nuageuses et des contenus en eau
    29772988
    2978      if (iflag_cldcon>=5) then
     2989     if (iflag_cldth>=5) then
    29792990
    29802991        do k=1,klev
    … …  
    31323143     calday = REAL(days_elapsed + 1) + jH_cur
    31333144
    3134      call chemtime(itap+itau_phy-1, date0, dtime)
     3145     call chemtime(itap+itau_phy-1, date0, dtime, itap)
    31353146     IF (config_inca == 'aero' .OR. config_inca == 'aeNP') THEN
    31363147        CALL AEROSOL_METEO_CALC( &
    … …  
    34553466  ! Ajouter la tendance des rayonnements (tous les pas)
    34563467  !
    3457   DO k = 1, klev
    3458      DO i = 1, klon
    3459         t_seri(i,k) = t_seri(i,k) &
    3460              + (heat(i,k)-cool(i,k)) * dtime/RDAY
    3461      ENDDO
    3462   ENDDO
     3468  d_t_swr(:,:)=heat(:,:)*dtime/RDAY
     3469  d_t_lwr(:,:)=-cool(:,:)*dtime/RDAY
     3470  d_t_sw0(:,:)=heat0(:,:)*dtime/RDAY
     3471  d_t_lw0(:,:)=-cool0(:,:)*dtime/RDAY
     3472  CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW')
     3473  CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW')
     3474
    34633475  !
    34643476  if (mydebug) then

  • LMDZ5/branches/testing/libf/phylmd/phytrac_mod.F90

    r2187 r2220  
    310310    !RomP >>>
    311311    INTEGER,SAVE  :: iflag_lscav_omp,iflag_lscav
     312    REAL, SAVE ::   ccntrAA_in,ccntrAA_omp
     313    REAL, SAVE ::   ccntrENV_in,ccntrENV_omp
     314    REAL, SAVE ::   coefcoli_in,coefcoli_omp
     315
    312316    LOGICAL,SAVE  :: convscav_omp,convscav
    313317!$OMP THREADPRIVATE(iflag_lscav)
     318!$OMP THREADPRIVATE(ccntrAA_in,ccntrENV_in,coefcoli_in)
    314319!$OMP THREADPRIVATE(convscav)
    315320    !RomP <<<
    … …  
    412417       iflag_lscav_omp=1
    413418       call getin('iflag_lscav', iflag_lscav_omp)
     419       ccntrAA_omp=1
     420       ccntrENV_omp=1.
     421       coefcoli_omp=0.001
     422       call getin('ccntrAA', ccntrAA_omp)
     423       call getin('ccntrENV', ccntrENV_omp)
     424       call getin('coefcoli', coefcoli_omp)
    414425!$OMP END MASTER
    415426!$OMP BARRIER
    416427       iflag_lscav=iflag_lscav_omp
     428       ccntrAA_in=ccntrAA_omp
     429       ccntrENV_in=ccntrENV_omp
     430       coefcoli_in=coefcoli_omp
    417431       !
    418432       SELECT CASE(iflag_lscav)
    … …  
    463477                IF (convscav.and.aerosol(it)) THEN
    464478                   flag_cvltr(it)=.true.
    465                    ccntrAA(it) =1.0         !--a modifier par JYG a lire depuis fichier
    466                    ccntrENV(it)=1.0
    467                    coefcoli(it)=0.001
     479                   ccntrAA(it) =ccntrAA_in    !--a modifier par JYG a lire depuis fichier
     480                   ccntrENV(it)=ccntrENV_in
     481                   coefcoli(it)=coefcoli_in
    468482                ELSE
    469483                   flag_cvltr(it)=.false.
    … …  
    613627                !--with the full array tr_seri even if only item it is processed
    614628
     629                print*,'CV SCAV ',it,ccntrAA(it),ccntrENV(it)
     630
    615631                CALL cvltr_scav(pdtphys, da, phi,phi2,d1a,dam, mp,ep,         &
    616632                     sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm,              &     
    … …  
    747763          !
    748764          DO it = 1, nbtr
     765
     766             IF (aerosol(it)) THEN
    749767             !  incloud scavenging and removal by large scale rain ! orig : ql_incl was replaced by 0.5e-3 kg/kg
    750768             ! the value 0.5e-3 kg/kg is from Giorgi and Chameides (1986), JGR
    … …  
    763781             ENDDO
    764782             CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'lsc scav it = '//solsym(it))
     783             ENDIF
     784
    765785          END DO  !tr
    766786

  • LMDZ5/branches/testing/libf/phylmd/radlwsw_m.F90

    r2160 r2220  
    765765! RII0 = RIP0M15 ! =rip0m if Morcrette non-each time step call.
    766766         RII0=solaire/zdist/zdist
    767         print*,'+++ radlwsw: solaire ,RII0',solaire,RII0
     767!print*,'+++ radlwsw: solaire ,RII0',solaire,RII0
    768768!  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    769769! Ancien appel a RECMWF (celui du cy25)
    … …  
    819819         ok_ade, ok_aie, flag_aerosol,flag_aerosol_strat) ! flags aerosols
    820820           
    821          print *,'RADLWSW: apres RECMWF'
     821!        print *,'RADLWSW: apres RECMWF'
    822822      if(lldebug) then
    823823        CALL writefield_phy('zemtd_i',ZEMTD_i,klev+1)
    … …  
    954954         zsollwdown(i)= ZFLDN(i,1)
    955955      ENDDO
    956       print*,'OK2'
     956!     print*,'OK2'
    957957
    958958! extrait de SW_AR4

  • LMDZ5/branches/testing/libf/phylmd/rrtm/radlsw.F90

    r2160 r2220  
    10621062!                 ------------------------------------
    10631063
    1064 print *,'RADLSW: LPHYLIN, LRRTM',LPHYLIN, LRRTM
     1064!print *,'RADLSW: LPHYLIN, LRRTM',LPHYLIN, LRRTM
    10651065IF (.NOT.LPHYLIN) THEN
    10661066  IF ( .NOT. LRRTM) THEN
    … …  
    10741074     & ZEMIT , PFLUX , PFLUC &
    10751075     & ) 
    1076     print *,'RADLSW: apres CALL LW'
     1076!   print *,'RADLSW: apres CALL LW'
    10771077    IF(LLDEBUG) THEN
    10781078    call writefield_phy('radlsw_flux1',PFLUX(:,1,:),klev+1)
    … …  
    11011101    ENDDO
    11021102
    1103     print *,'RADLSW: avant CALL RRTM_RRTM_140GP,PAP=',PAP(1,:)
     1103!   print *,'RADLSW: avant CALL RRTM_RRTM_140GP,PAP=',PAP(1,:)
    11041104    CALL RRTM_RRTM_140GP &
    11051105     & ( KIDIA , KFDIA , KLON  , KLEV,&
    … …  
    11111111     & PTAU_LW,&
    11121112     & ZEMIT , PFLUX , PFLUC , ZTCLEAR )
    1113     print *,'RADLSW: apres CALL RRTM_RRTM_140GP'
     1113!   print *,'RADLSW: apres CALL RRTM_RRTM_140GP'
    11141114
    11151115  ENDIF
    … …  
    11181118  PFLUX(:,:,:)= 0.0_JPRB
    11191119  PFLUC(:,:,:)= 0.0_JPRB
    1120   print *,'RADLSW: ZEMIT,PFLUX et PFLUC = 0'
     1120! print *,'RADLSW: ZEMIT,PFLUX et PFLUC = 0'
    11211121ENDIF
    11221122

  • LMDZ5/branches/testing/libf/phylmd/rrtm/sw1s.F90

    r1999 r2220  
    304304
    305305ELSEIF (NSW == 6) THEN
    306 print *,'... dans SW1S: NSW=',NSW
     306!print *,'... dans SW1S: NSW=',NSW
    307307
    308308!*         3.2   SIX SPECTRAL INTERVALS

  • LMDZ5/branches/testing/libf/phylmd/surf_land_mod.F90

    r1910 r2220  
    8585    REAL, DIMENSION(klon) :: pref_tmp
    8686    REAL, DIMENSION(klon) :: swdown     ! downwelling shortwave radiation at land surface
    87     REAL, DIMENSION(klon) :: lwdown     ! downwelling longwave radiation at land surface
    8887    REAL, DIMENSION(klon) :: epot_air           ! potential air temperature
    8988    REAL, DIMENSION(klon) :: tsol_rad, emis_new ! output from interfsol not used
    … …  
    106105       pref_tmp(1:knon)  = pref(1:knon)/100.
    107106!
    108 !* Calculate incoming flux for SW and LW interval: swdown, lwdown
     107!* Calculate incoming flux for SW and LW interval: swdown
    109108!
    110109       swdown(:) = 0.0
    111        lwdown(:) = 0.0
    112110       DO i = 1, knon
    113111          swdown(i) = swnet(i)/(1-albedo(i))
    114           lwdown(i) = lwnet(i) + RSIGMA*tsurf(i)**4
    115112       END DO
    116113!

  • LMDZ5/branches/testing/libf/phylmd/surf_ocean_mod.F90

    r2187 r2220  
    123123            AcoefU, AcoefV, BcoefU, BcoefV, &
    124124            ps, u1, v1, tsurf_in, &
    125             radsol, snow, agesno, &
     125            radsol, snow, &
    126126            qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
    127127            tsurf_new, dflux_s, dflux_l, lmt_bils)

  • LMDZ5/branches/testing/libf/phylmd/surf_seaice_mod.F90

    r2073 r2220  
    2525  USE ocean_forced_mod, ONLY : ocean_forced_ice
    2626  USE ocean_cpl_mod, ONLY    : ocean_cpl_ice
     27  USE ocean_slab_mod, ONLY   : ocean_slab_ice
    2728  USE indice_sol_mod
    2829
    … …  
    108109            tsurf_new, dflux_s, dflux_l)
    109110       
    110 !    ELSE IF (type_ocean == 'slab'.AND.version_ocean=='sicINT') THEN
    111 !       CALL ocean_slab_ice( &
    112 !          itime, dtime, jour, knon, knindex, &
    113 !          debut, &
    114 !          tsurf, p1lay, cdragh, precip_rain, precip_snow, temp_air, spechum,&
    115 !          AcoefH, AcoefQ, BcoefH, BcoefQ, &
    116 !          ps, u1, v1, &
    117 !          radsol, snow, qsurf, qsol, agesno, tsoil, &
    118 !          alb1_new, alb2_new, evap, fluxsens, fluxlat, &
    119 !          tsurf_new, dflux_s, dflux_l)
    120 !
     111    ELSE IF (type_ocean == 'slab'.AND.version_ocean=='sicINT') THEN
     112       CALL ocean_slab_ice( &
     113          itime, dtime, jour, knon, knindex, &
     114          tsurf, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum,&
     115          AcoefH, AcoefQ, BcoefH, BcoefQ, &
     116            AcoefU, AcoefV, BcoefU, BcoefV, &
     117          ps, u1, v1, &
     118          radsol, snow, qsurf, qsol, agesno, &
     119          alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
     120          tsurf_new, dflux_s, dflux_l, swnet)
     121
    121122      ELSE ! type_ocean=force or slab +sicOBS or sicNO
    122123       CALL ocean_forced_ice( &

  • LMDZ5/branches/testing/libf/phylmd/surface_data.F90

    r2160 r2220  
    77  REAL, PARAMETER        :: tau_gl=86400.*5.
    88  REAL, PARAMETER        :: calsno=1./(2.3867e+06*.15)
    9   REAL, PARAMETER        :: t_freeze=271.35
    10   REAL, PARAMETER        :: t_melt=273.15
    119 
    1210  LOGICAL, SAVE          :: ok_veget      ! true for use of vegetation model ORCHIDEE

  • LMDZ5/branches/testing/libf/phylmd/thermcell_main.F90

    r2160 r2220  
    976976      enddo
    977977
     978! Ale sec (max de wmax/2 sous la zone d'inhibition) dans
     979! le cas iflag_trig_bl=3
     980      IF (iflag_trig_bl==3) Ale_bl(:)=0.5*wmax_sec(:)**2
     981
    978982!test:calcul de la ponderation des couches pour KE
    979983!initialisations

  • LMDZ5/branches/testing/libf/phylmd/thermcell_plume.F90

    r2187 r2220  
    359359                   zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
    360360    &                    +(1.-fact_shell)*zbuoy(ig,l)
    361          print*,'on est pass?? par l??1',l,lt,ztv1,ztv2,ztv(ig,lt),ztv_est(ig,l),ztva_est(ig,l),ztv(ig,l), &
    362    &                                   zinv,zlmelup,zbuoy(ig,l),zbuoyjam(ig,l)
    363361                elseif (zlmelup.ge.zinv) then
    364362                 ztv_est2=atv2*0.5*(zlmelup+zinv)+btv2
    … …  
    370368    &            ztv_est1)/ztv_est1)+(1.-fact_shell)*zbuoy(ig,l)
    371369
    372          print*,'on est pass?? par l??2',l,lt,ztv_est1,ztv_est2,ztv(ig,lt),ztv_est(ig,l),ztva_est(ig,l),ztv(ig,l), &
    373    &                                   zinv,zlmelup,zbuoy(ig,l),zbuoyjam(ig,l)
    374370                else
    375371                   ztv_est(ig,l)=atv1*zlmel+btv1
    376372                   zbuoyjam(ig,l)=fact_shell*RG*(ztva_est(ig,l)-ztv_est(ig,l))/ztv_est(ig,l) &
    377373    &                           +(1.-fact_shell)*zbuoy(ig,l)
    378                  print*,'on est pass?? par l??3',l,lt,ztv1,ztv2,ztv(ig,lt),ztv_est(ig,l),ztva_est(ig,l),ztv(ig,l), &
    379    &                                   zinv,zlmelup,zbuoy(ig,l),zbuoyjam(ig,l)
    380374                endif
    381375
    … …  
    392386                endif
    393387
    394                  print*,'on est pass?? par l??4',l,lt,ztv1,ztv2,ztv(ig,lt),ztv(ig,l),ztva_est(ig,l), &
    395    &                                   zlmelup,zbuoy(ig,l),zbuoyjam(ig,l)
    396388!          zbuoyjam(ig,l)=fact_shell*RG*(((zlmelup-zltdwn)/zdz)*(ztva_est(ig,l)- &
    397389!    &          ztv1)/ztv1+((zltdwn-zlmeldwn)/zdz)*(ztva_est(ig,l)- &
    … …  
    417409    &          ztv(ig,lt))/ztv(ig,lt)+((zdz2-lmel)/zdz3)*(ztva_est(ig,l)- &
    418410    &          ztv(ig,lt-1))/ztv(ig,lt-1))+0.*zbuoy(ig,l)
    419          print*,'on est pass?? par l??',l,lt,ztv(ig,lt),ztva_est(ig,l),ztv(ig,l), &
    420    &                                   zbuoy(ig,l),zbuoyjam(ig,l)
    421411        endif !   if (iflag_thermals_ed.lt.8) then
    422412

  • LMDZ5/branches/testing/libf/phylmd/tilft43.F90

    r1999 r2220  
    33
    44SUBROUTINE tlift43(p, t, q, qs, gz, icb, nk, tvp, tpk, clw, nd, nl, kk)
     5  IMPLICIT NONE
    56  REAL gz(nd), tpk(nd), clw(nd), p(nd)
    67  REAL t(nd), q(nd), qs(nd), tvp(nd), lv0
    7 
     8  INTEGER icb, nk, nd, nl, kk
     9  REAL cpd, cpv,  cl, g, rowl, gravity, cpvmcl, eps, epsi
     10  REAL ah0, cpp, cpinv, tg, qg, alv, s, ahg, tc, denom, es
     11  INTEGER i, nst, nsb, j
    812  ! ***   ASSIGN VALUES OF THERMODYNAMIC CONSTANTS     ***
    913

  • LMDZ5/branches/testing/libf/phylmd/tlift.F90

    r1999 r2220  
    44SUBROUTINE tlift(p, t, rr, rs, gz, plcl, icb, nk, tvp, tpk, clw, nd, nl, &
    55    dtvpdt1, dtvpdq1)
    6 
     6  IMPLICIT NONE
    77  ! Argument NK ajoute (jyg) = Niveau de depart de la
    88  ! convection
    9 
    10   PARAMETER (na=60)
    11   REAL gz(nd), tpk(nd), clw(nd)
     9  INTEGER icb, nk, nd, nl
     10  INTEGER,PARAMETER :: na=60
     11  REAL gz(nd), tpk(nd), clw(nd), plcl
    1212  REAL t(nd), rr(nd), rs(nd), tvp(nd), p(nd)
    1313  REAL dtvpdt1(nd), dtvpdq1(nd) ! Derivatives of parcel virtual
    … …  
    1717  REAL dtpdt1(na), dtpdq1(na) ! Derivatives of parcel temperature
    1818  ! wrt T1 and Q1
    19 
     19  REAL gravity, cpd, cpv, cl, ci, cpvmcl, clmci, eps, alv0, alf0
     20  REAL cpp, cpinv, ah0, alf, tg, s, ahg, tc, denom, alv, es, esi
     21  REAL qsat_new, snew
     22  INTEGER icbl, i, imin, j, icb1
    2023
    2124  LOGICAL ice_conv

  • LMDZ5/branches/testing/libf/phylmd/wake.F90

    r2160 r2220  
    17561756  ! a une humidite positive dans la zone (x) et dans la zone (w).
    17571757  ! ------------------------------------------------------
    1758 
     1758  IMPLICIT NONE
    17591759
    17601760  ! Input
    … …  
    17721772  REAL epsilon
    17731773  ! DATA epsilon/1.e-15/
     1774  INTEGER i,k
    17741775
    17751776  DO k = 1, nl
Note: See TracChangeset for help on using the changeset viewer.