Ignore:
Timestamp:
Dec 23, 2010, 5:38:42 PM (14 years ago)
Author:
lguez
Message:

Conversion to free source form for "disvert" and "iniacademic", no
other change. Bug fix in "fisrtilp": "fraca" could appear in an
expression while not defined.

Location:
LMDZ5/branches/LMDZ5V2.0-dev/libf
Files:
5 moved

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/disvert.F90

    r1463 r1472  
    1 !
    21! $Id$
    3 !
    4       SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
    52
    6 c    Auteur :  P. Le Van .
    7 c
    8       IMPLICIT NONE
     3SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
    94
    10 #include "dimensions.h"
    11 #include "paramet.h"
    12 #include "iniprint.h"
    13 #include "logic.h"
    14 c
    15 c=======================================================================
    16 c
    17 c
    18 c    s = sigma ** kappa   :  coordonnee  verticale
    19 c    dsig(l)            : epaisseur de la couche l ds la coord.  s
    20 c    sig(l)             : sigma a l'interface des couches l et l-1
    21 c    ds(l)              : distance entre les couches l et l-1 en coord.s
    22 c
    23 c=======================================================================
    24 c
    25       REAL pa,preff
    26       REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
    27       REAL presnivs(llm)
    28 c
    29 c   declarations:
    30 c   -------------
    31 c
    32       REAL sig(llm+1),dsig(llm)
    33        real zzz(1:llm+1)
    34        real dzz(1:llm)
    35       real zk,zkm1,dzk1,dzk2,k0,k1
    36 c
    37       INTEGER l
    38       REAL snorm,dsigmin
    39       REAL alpha,beta,gama,delta,deltaz,h
    40       INTEGER np,ierr
    41       REAL pi,x
     5  ! Auteur : P. Le Van
    426
    43       REAL SSUM
    44 c
    45 c-----------------------------------------------------------------------
    46 c
    47       pi=2.*ASIN(1.)
     7  IMPLICIT NONE
    488
    49       OPEN(99,file='sigma.def',status='old',form='formatted',
    50      s   iostat=ierr)
     9  include "dimensions.h"
     10  include "paramet.h"
     11  include "iniprint.h"
     12  include "logic.h"
    5113
    52 c-----------------------------------------------------------------------
    53 c   cas 1 on lit les options dans sigma.def:
    54 c   ----------------------------------------
     14  ! s = sigma ** kappa : coordonnee verticale
     15  ! dsig(l) : epaisseur de la couche l ds la coord. s
     16  ! sig(l) : sigma a l'interface des couches l et l-1
     17  ! ds(l) : distance entre les couches l et l-1 en coord.s
    5518
    56       IF (ierr.eq.0) THEN
     19  REAL pa, preff
     20  REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1)
     21  REAL presnivs(llm)
    5722
    58       READ(99,*) h           ! hauteur d'echelle 8.
    59       READ(99,*) deltaz      ! epaiseur de la premiere couche 0.04
    60       READ(99,*) beta        ! facteur d'acroissement en haut 1.3
    61       READ(99,*) k0          ! nombre de couches dans la transition surf
    62       READ(99,*) k1          ! nombre de couches dans la transition haute
    63       CLOSE(99)
    64       alpha=deltaz/(llm*h)
    65       write(lunout,*)'h,alpha,k0,k1,beta'
     23  REAL sig(llm+1), dsig(llm)
     24  real zk, zkm1, dzk1, dzk2, k0, k1
    6625
    67 c     read(*,*) h,deltaz,beta,k0,k1 ! 8 0.04 4 20 1.2
     26  INTEGER l
     27  REAL dsigmin
     28  REAL alpha, beta, deltaz, h
     29  INTEGER iostat
     30  REAL pi, x
    6831
    69       alpha=deltaz/tanh(1./k0)*2.
    70       zkm1=0.
    71       sig(1)=1.
    72       do l=1,llm
    73         sig(l+1)=(cosh(l/k0))**(-alpha*k0/h)
    74      + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
     32  !-----------------------------------------------------------------------
     33
     34  pi = 2 * ASIN(1.)
     35
     36  OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
     37
     38  IF (iostat == 0) THEN
     39     ! cas 1 on lit les options dans sigma.def:
     40     READ(99, *) h ! hauteur d'echelle 8.
     41     READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
     42     READ(99, *) beta ! facteur d'acroissement en haut 1.3
     43     READ(99, *) k0 ! nombre de couches dans la transition surf
     44     READ(99, *) k1 ! nombre de couches dans la transition haute
     45     CLOSE(99)
     46     alpha=deltaz/(llm*h)
     47     write(lunout, *)'h, alpha, k0, k1, beta'
     48
     49     alpha=deltaz/tanh(1./k0)*2.
     50     zkm1=0.
     51     sig(1)=1.
     52     do l=1, llm
     53        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) &
     54             *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
    7555        zk=-h*log(sig(l+1))
    7656
    7757        dzk1=alpha*tanh(l/k0)
    7858        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
    79         write(lunout,*)l,sig(l+1),zk,zk-zkm1,dzk1,dzk2
     59        write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
    8060        zkm1=zk
    81       enddo
     61     enddo
    8262
    83       sig(llm+1)=0.
     63     sig(llm+1)=0.
    8464
    85 c
    86        DO 2  l = 1, llm
    87        dsig(l) = sig(l)-sig(l+1)
    88    2   CONTINUE
    89 c
     65     DO l = 1, llm
     66        dsig(l) = sig(l)-sig(l+1)
     67     end DO
     68  ELSE
     69     if (ok_strato) then
     70        if (llm==39) then
     71           dsigmin=0.3
     72        else if (llm==50) then
     73           dsigmin=1.
     74        else
     75           WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
     76           dsigmin=1.
     77        endif
     78        WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
     79     endif
    9080
    91       ELSE
    92 c-----------------------------------------------------------------------
    93 c   cas 2 ancienne discretisation (LMD5...):
    94 c   ----------------------------------------
     81     DO l = 1, llm
     82        x = pi * (l - 0.5) / (llm + 1)
    9583
    96       WRITE(LUNOUT,*)'WARNING!!! Ancienne discretisation verticale'
     84        IF (ok_strato) THEN
     85           dsig(l) =(dsigmin + 7. * SIN(x)**2) &
     86                *(1. - tanh(2 * x / pi - 1.))**2 / 4.
     87        ELSE
     88           dsig(l) = 1.0 + 7.0 * SIN(x)**2
     89        ENDIF
     90     ENDDO
     91     dsig = dsig / sum(dsig)
     92     sig(llm+1) = 0.
     93     DO l = llm, 1, -1
     94        sig(l) = sig(l+1) + dsig(l)
     95     ENDDO
     96  ENDIF
    9797
    98       if (ok_strato) then
    99          if (llm==39) then
    100             dsigmin=0.3
    101          else if (llm==50) then
    102             dsigmin=1.
    103          else
    104             WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
    105             dsigmin=1.
    106          endif
    107          WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
    108       endif
     98  DO l=1, llm
     99     nivsigs(l) = REAL(l)
     100  ENDDO
    109101
    110       h=7.
    111       snorm  = 0.
    112       DO l = 1, llm
    113          x = 2.*asin(1.) * (REAL(l)-0.5) / REAL(llm+1)
     102  DO l=1, llmp1
     103     nivsig(l)= REAL(l)
     104  ENDDO
    114105
    115          IF (ok_strato) THEN
    116            dsig(l) =(dsigmin + 7.0 * SIN(x)**2)
    117      &            *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2       
    118          ELSE
    119            dsig(l) = 1.0 + 7.0 * SIN(x)**2
    120          ENDIF
     106  ! .... Calculs de ap(l) et de bp(l) ....
     107  ! ..... pa et preff sont lus sur les fichiers start par lectba .....
    121108
    122          snorm = snorm + dsig(l)
    123       ENDDO
    124       snorm = 1./snorm
    125       DO l = 1, llm
    126          dsig(l) = dsig(l)*snorm
    127       ENDDO
    128       sig(llm+1) = 0.
    129       DO l = llm, 1, -1
    130          sig(l) = sig(l+1) + dsig(l)
    131       ENDDO
     109  bp(llmp1) = 0.
    132110
    133       ENDIF
     111  DO l = 1, llm
     112     bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
     113     ap(l) = pa * ( sig(l) - bp(l) )
     114  ENDDO
    134115
     116  bp(1)=1.
     117  ap(1)=0.
    135118
    136       DO l=1,llm
    137         nivsigs(l) = REAL(l)
    138       ENDDO
     119  ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
    139120
    140       DO l=1,llmp1
    141         nivsig(l)= REAL(l)
    142       ENDDO
     121  write(lunout, *)' BP '
     122  write(lunout, *) bp
     123  write(lunout, *)' AP '
     124  write(lunout, *) ap
    143125
    144 c
    145 c    ....  Calculs  de ap(l) et de bp(l)  ....
    146 c    .........................................
    147 c
    148 c
    149 c   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
    150 c
     126  write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
     127  write(lunout, *)'couches calcules pour une pression de surface =', preff
     128  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de'
     129  write(lunout, *)'8km'
     130  DO l = 1, llm
     131     dpres(l) = bp(l) - bp(l+1)
     132     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
     133     write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
     134          log(preff/presnivs(l))*8. &
     135          , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ &
     136          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
     137  ENDDO
    151138
    152       bp(llmp1) =   0.
     139  write(lunout, *) 'PRESNIVS '
     140  write(lunout, *) presnivs
    153141
    154       DO l = 1, llm
    155 cc
    156 ccc    ap(l) = 0.
    157 ccc    bp(l) = sig(l)
    158 
    159       bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
    160       ap(l) = pa * ( sig(l) - bp(l) )
    161 c
    162       ENDDO
    163 
    164       bp(1)=1.
    165       ap(1)=0.
    166 
    167       ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
    168 
    169       write(lunout,*)' BP '
    170       write(lunout,*)  bp
    171       write(lunout,*)' AP '
    172       write(lunout,*)  ap
    173 
    174       write(lunout,*)
    175      .'Niveaux de pressions approximatifs aux centres des'
    176       write(lunout,*)'couches calcules pour une pression de surface =',
    177      .                 preff
    178       write(lunout,*)
    179      .     'et altitudes equivalentes pour une hauteur d echelle de'
    180       write(lunout,*)'8km'
    181       DO l = 1, llm
    182        dpres(l) = bp(l) - bp(l+1)
    183        presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
    184        write(lunout,*)'PRESNIVS(',l,')=',presnivs(l),'    Z ~ ',
    185      .        log(preff/presnivs(l))*8.
    186      .  ,'   DZ ~ ',8.*log((ap(l)+bp(l)*preff)/
    187      .       max(ap(l+1)+bp(l+1)*preff,1.e-10))
    188       ENDDO
    189 
    190       write(lunout,*)' PRESNIVS '
    191       write(lunout,*)presnivs
    192 
    193       RETURN
    194       END
     142END SUBROUTINE disvert

  • LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/iniacademic.F90

    r1458 r1472  
    22! $Id$
    33!
    4 c
    5 c
    6       SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
    7 
    8       USE filtreg_mod
    9       USE infotrac, ONLY : nqtot
    10       USE control_mod, ONLY: day_step,planet_type
     4SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
     5
     6  USE filtreg_mod
     7  USE infotrac, ONLY : nqtot
     8  USE control_mod, ONLY: day_step,planet_type
    119#ifdef CPP_IOIPSL
    12       USE IOIPSL
     10  USE IOIPSL
    1311#else
    14 ! if not using IOIPSL, we still need to use (a local version of) getin
    15       USE ioipsl_getincom
     12  ! if not using IOIPSL, we still need to use (a local version of) getin
     13  USE ioipsl_getincom
    1614#endif
    17       USE Write_Field
    18  
    19 
    20 c%W%    %G%
    21 c=======================================================================
    22 c
    23 c   Author:    Frederic Hourdin      original: 15/01/93
    24 c   -------
    25 c
    26 c   Subject:
    27 c   ------
    28 c
    29 c   Method:
    30 c   --------
    31 c
    32 c   Interface:
    33 c   ----------
    34 c
    35 c      Input:
    36 c      ------
    37 c
    38 c      Output:
    39 c      -------
    40 c
    41 c=======================================================================
    42       IMPLICIT NONE
    43 c-----------------------------------------------------------------------
    44 c   Declararations:
    45 c   ---------------
    46 
    47 #include "dimensions.h"
    48 #include "paramet.h"
    49 #include "comvert.h"
    50 #include "comconst.h"
    51 #include "comgeom.h"
    52 #include "academic.h"
    53 #include "ener.h"
    54 #include "temps.h"
    55 #include "iniprint.h"
    56 #include "logic.h"
    57 
    58 c   Arguments:
    59 c   ----------
    60 
    61       real time_0
    62 
    63 c   variables dynamiques
    64       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
    65       REAL teta(ip1jmp1,llm)                 ! temperature potentielle
    66       REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
    67       REAL ps(ip1jmp1)                       ! pression  au sol
    68       REAL masse(ip1jmp1,llm)                ! masse d'air
    69       REAL phis(ip1jmp1)                     ! geopotentiel au sol
    70 
    71 c   Local:
    72 c   ------
    73 
    74       REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
    75       REAL pks(ip1jmp1)                      ! exner au  sol
    76       REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
    77       REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
    78       REAL phi(ip1jmp1,llm)                  ! geopotentiel
    79       REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
    80       real tetajl(jjp1,llm)
    81       INTEGER i,j,l,lsup,ij
    82 
    83       REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
    84       REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
    85       LOGICAL ok_geost             ! Initialisation vent geost. ou nul
    86       LOGICAL ok_pv                ! Polar Vortex
    87       REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex
    88      
    89       real zz,ran1
    90       integer idum
    91 
    92       REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
    93 
    94 c-----------------------------------------------------------------------
    95 ! 1. Initializations for Earth-like case
    96 ! --------------------------------------
    97 c
    98         ! initialize planet radius, rotation rate,...
    99         call conf_planete
    100 
    101         time_0=0.
    102         day_ref=1
    103         annee_ref=0
    104 
    105         im         = iim
    106         jm         = jjm
    107         day_ini    = 1
    108         dtvr    = daysec/REAL(day_step)
    109         zdtvr=dtvr
    110         etot0      = 0.
    111         ptot0      = 0.
    112         ztot0      = 0.
    113         stot0      = 0.
    114         ang0       = 0.
    115 
    116         if (llm.eq.1) then
    117           ! specific initializations for the shallow water case
    118           kappa=1
     15  USE Write_Field
     16
     17  !   Author:    Frederic Hourdin      original: 15/01/93
     18
     19  IMPLICIT NONE
     20
     21  !   Declararations:
     22  !   ---------------
     23
     24  include "dimensions.h"
     25  include "paramet.h"
     26  include "comvert.h"
     27  include "comconst.h"
     28  include "comgeom.h"
     29  include "academic.h"
     30  include "ener.h"
     31  include "temps.h"
     32  include "iniprint.h"
     33  include "logic.h"
     34
     35  !   Arguments:
     36  !   ----------
     37
     38  real time_0
     39
     40  !   variables dynamiques
     41  REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
     42  REAL teta(ip1jmp1,llm)                 ! temperature potentielle
     43  REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
     44  REAL ps(ip1jmp1)                       ! pression  au sol
     45  REAL masse(ip1jmp1,llm)                ! masse d'air
     46  REAL phis(ip1jmp1)                     ! geopotentiel au sol
     47
     48  !   Local:
     49  !   ------
     50
     51  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
     52  REAL pks(ip1jmp1)                      ! exner au  sol
     53  REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
     54  REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
     55  REAL phi(ip1jmp1,llm)                  ! geopotentiel
     56  REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
     57  real tetajl(jjp1,llm)
     58  INTEGER i,j,l,lsup,ij
     59
     60  REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
     61  REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
     62  LOGICAL ok_geost             ! Initialisation vent geost. ou nul
     63  LOGICAL ok_pv                ! Polar Vortex
     64  REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex
     65
     66  real zz,ran1
     67  integer idum
     68
     69  REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
     70
     71  !-----------------------------------------------------------------------
     72  ! 1. Initializations for Earth-like case
     73  ! --------------------------------------
     74  !
     75  ! initialize planet radius, rotation rate,...
     76  call conf_planete
     77
     78  time_0=0.
     79  day_ref=1
     80  annee_ref=0
     81
     82  im         = iim
     83  jm         = jjm
     84  day_ini    = 1
     85  dtvr    = daysec/REAL(day_step)
     86  zdtvr=dtvr
     87  etot0      = 0.
     88  ptot0      = 0.
     89  ztot0      = 0.
     90  stot0      = 0.
     91  ang0       = 0.
     92
     93  if (llm == 1) then
     94     ! specific initializations for the shallow water case
     95     kappa=1
     96  endif
     97
     98  CALL iniconst
     99  CALL inigeom
     100  CALL inifilr
     101
     102  if (llm == 1) then
     103     ! initialize fields for the shallow water case, if required
     104     if (.not.read_start) then
     105        phis(:)=0.
     106        q(:,:,:)=0
     107        CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
     108     endif
     109  endif
     110
     111  academic_case: if (iflag_phys == 2) then
     112     ! initializations
     113
     114     ! 1. local parameters
     115     ! by convention, winter is in the southern hemisphere
     116     ! Geostrophic wind or no wind?
     117     ok_geost=.TRUE.
     118     CALL getin('ok_geost',ok_geost)
     119     ! Constants for Newtonian relaxation and friction
     120     k_f=1.                !friction
     121     CALL getin('k_j',k_f)
     122     k_f=1./(daysec*k_f)
     123     k_c_s=4.  !cooling surface
     124     CALL getin('k_c_s',k_c_s)
     125     k_c_s=1./(daysec*k_c_s)
     126     k_c_a=40. !cooling free atm
     127     CALL getin('k_c_a',k_c_a)
     128     k_c_a=1./(daysec*k_c_a)
     129     ! Constants for Teta equilibrium profile
     130     teta0=315.     ! mean Teta (S.H. 315K)
     131     CALL getin('teta0',teta0)
     132     ttp=200.       ! Tropopause temperature (S.H. 200K)
     133     CALL getin('ttp',ttp)
     134     eps=0.         ! Deviation to N-S symmetry(~0-20K)
     135     CALL getin('eps',eps)
     136     delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
     137     CALL getin('delt_y',delt_y)
     138     delt_z=10.     ! Vertical Gradient (S.H. 10K)
     139     CALL getin('delt_z',delt_z)
     140     ! Polar vortex
     141     ok_pv=.false.
     142     CALL getin('ok_pv',ok_pv)
     143     phi_pv=-50.            ! Latitude of edge of vortex
     144     CALL getin('phi_pv',phi_pv)
     145     phi_pv=phi_pv*pi/180.
     146     dphi_pv=5.             ! Width of the edge
     147     CALL getin('dphi_pv',dphi_pv)
     148     dphi_pv=dphi_pv*pi/180.
     149     gam_pv=4.              ! -dT/dz vortex (in K/km)
     150     CALL getin('gam_pv',gam_pv)
     151
     152     ! 2. Initialize fields towards which to relax
     153     ! Friction
     154     knewt_g=k_c_a
     155     DO l=1,llm
     156        zsig=presnivs(l)/preff
     157        knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
     158        kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
     159     ENDDO
     160     DO j=1,jjp1
     161        clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
     162     ENDDO
     163
     164     ! Potential temperature
     165     DO l=1,llm
     166        zsig=presnivs(l)/preff
     167        tetastrat=ttp*zsig**(-kappa)
     168        tetapv=tetastrat
     169        IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
     170           tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
     171        ENDIF
     172        DO j=1,jjp1
     173           ! Troposphere
     174           ddsin=sin(rlatu(j))
     175           tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
     176                -delt_z*(1.-ddsin*ddsin)*log(zsig)
     177           ! Profil stratospherique isotherme (+vortex)
     178           w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
     179           tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
     180           tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
     181        ENDDO
     182     ENDDO
     183
     184     !          CALL writefield('theta_eq',tetajl)
     185
     186     do l=1,llm
     187        do j=1,jjp1
     188           do i=1,iip1
     189              ij=(j-1)*iip1+i
     190              tetarappel(ij,l)=tetajl(j,l)
     191           enddo
     192        enddo
     193     enddo
     194
     195     ! 3. Initialize fields (if necessary)
     196     IF (.NOT. read_start) THEN
     197        ! surface pressure
     198        ps(:)=preff
     199        ! ground geopotential
     200        phis(:)=0.
     201
     202        CALL pression ( ip1jmp1, ap, bp, ps, p       )
     203        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     204        CALL massdair(p,masse)
     205
     206        ! bulk initialization of temperature
     207        teta(:,:)=tetarappel(:,:)
     208
     209        ! geopotential
     210        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
     211
     212        ! winds
     213        if (ok_geost) then
     214           call ugeostr(phi,ucov)
     215        else
     216           ucov(:,:)=0.
    119217        endif
    120        
    121         CALL iniconst
    122         CALL inigeom
    123         CALL inifilr
    124 
    125         if (llm.eq.1) then
    126           ! initialize fields for the shallow water case, if required
    127           if (.not.read_start) then
    128             phis(:)=0.
    129             q(:,:,:)=0
    130             CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
    131           endif
    132         endif
    133 
    134         if (iflag_phys.eq.2) then
    135           ! initializations for the academic case
    136          
    137 !         if (planet_type=="earth") then
    138 
    139           ! 1. local parameters
    140           ! by convention, winter is in the southern hemisphere
    141           ! Geostrophic wind or no wind?
    142           ok_geost=.TRUE.
    143           CALL getin('ok_geost',ok_geost)
    144           ! Constants for Newtonian relaxation and friction
    145           k_f=1.                !friction
    146           CALL getin('k_j',k_f)
    147           k_f=1./(daysec*k_f)
    148           k_c_s=4.  !cooling surface
    149           CALL getin('k_c_s',k_c_s)
    150           k_c_s=1./(daysec*k_c_s)
    151           k_c_a=40. !cooling free atm
    152           CALL getin('k_c_a',k_c_a)
    153           k_c_a=1./(daysec*k_c_a)
    154           ! Constants for Teta equilibrium profile
    155           teta0=315.     ! mean Teta (S.H. 315K)
    156           CALL getin('teta0',teta0)
    157           ttp=200.       ! Tropopause temperature (S.H. 200K)
    158           CALL getin('ttp',ttp)
    159           eps=0.         ! Deviation to N-S symmetry(~0-20K)
    160           CALL getin('eps',eps)
    161           delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
    162           CALL getin('delt_y',delt_y)
    163           delt_z=10.     ! Vertical Gradient (S.H. 10K)
    164           CALL getin('delt_z',delt_z)
    165           ! Polar vortex
    166           ok_pv=.false.
    167           CALL getin('ok_pv',ok_pv)
    168           phi_pv=-50.            ! Latitude of edge of vortex
    169           CALL getin('phi_pv',phi_pv)
    170           phi_pv=phi_pv*pi/180.
    171           dphi_pv=5.             ! Width of the edge
    172           CALL getin('dphi_pv',dphi_pv)
    173           dphi_pv=dphi_pv*pi/180.
    174           gam_pv=4.              ! -dT/dz vortex (in K/km)
    175           CALL getin('gam_pv',gam_pv)
    176          
    177           ! 2. Initialize fields towards which to relax
    178           ! Friction
    179           knewt_g=k_c_a
    180           DO l=1,llm
    181            zsig=presnivs(l)/preff
    182            knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
    183            kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
    184           ENDDO
    185           DO j=1,jjp1
    186             clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
    187           ENDDO
    188          
    189           ! Potential temperature
    190           DO l=1,llm
    191            zsig=presnivs(l)/preff
    192            tetastrat=ttp*zsig**(-kappa)
    193            tetapv=tetastrat
    194            IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
    195                tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
    196            ENDIF
    197            DO j=1,jjp1
    198              ! Troposphere
    199              ddsin=sin(rlatu(j))
    200              tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin
    201      &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
    202              ! Profil stratospherique isotherme (+vortex)
    203              w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
    204              tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
    205              tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
    206            ENDDO
    207           ENDDO ! of DO l=1,llm
    208  
    209 !          CALL writefield('theta_eq',tetajl)
    210 
    211           do l=1,llm
    212             do j=1,jjp1
    213               do i=1,iip1
    214                  ij=(j-1)*iip1+i
    215                  tetarappel(ij,l)=tetajl(j,l)
    216               enddo
    217             enddo
    218           enddo
    219 
    220 
    221 !         else
    222 !          write(lunout,*)"iniacademic: planet types other than earth",
    223 !     &                   " not implemented (yet)."
    224 !          stop
    225 !         endif ! of if (planet_type=="earth")
    226 
    227           ! 3. Initialize fields (if necessary)
    228           IF (.NOT. read_start) THEN
    229             ! surface pressure
    230             ps(:)=preff
    231             ! ground geopotential
    232             phis(:)=0.
    233            
    234             CALL pression ( ip1jmp1, ap, bp, ps, p       )
    235             CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    236             CALL massdair(p,masse)
    237 
    238             ! bulk initialization of temperature
    239             teta(:,:)=tetarappel(:,:)
    240            
    241             ! geopotential
    242             CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
    243            
    244             ! winds
    245             if (ok_geost) then
    246               call ugeostr(phi,ucov)
    247             else
    248               ucov(:,:)=0.
    249             endif
    250             vcov(:,:)=0.
    251            
    252             ! bulk initialization of tracers
    253             if (planet_type=="earth") then
    254               ! Earth: first two tracers will be water
    255               do i=1,nqtot
    256                 if (i.eq.1) q(:,:,i)=1.e-10
    257                 if (i.eq.2) q(:,:,i)=1.e-15
    258                 if (i.gt.2) q(:,:,i)=0.
    259               enddo
    260             else
    261               q(:,:,:)=0
    262             endif ! of if (planet_type=="earth")
    263 
    264             ! add random perturbation to temperature
    265             idum  = -1
    266             zz = ran1(idum)
    267             idum  = 0
    268             do l=1,llm
    269               do ij=iip2,ip1jm
    270                 teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
    271               enddo
    272             enddo
    273 
    274             ! maintain periodicity in longitude
    275             do l=1,llm
    276               do ij=1,ip1jmp1,iip1
    277                 teta(ij+iim,l)=teta(ij,l)
    278               enddo
    279             enddo
    280 
    281           ENDIF ! of IF (.NOT. read_start)
    282         endif ! of if (iflag_phys.eq.2)
    283        
    284       END
    285 c-----------------------------------------------------------------------
     218        vcov(:,:)=0.
     219
     220        ! bulk initialization of tracers
     221        if (planet_type=="earth") then
     222           ! Earth: first two tracers will be water
     223           do i=1,nqtot
     224              if (i == 1) q(:,:,i)=1.e-10
     225              if (i == 2) q(:,:,i)=1.e-15
     226              if (i.gt.2) q(:,:,i)=0.
     227           enddo
     228        else
     229           q(:,:,:)=0
     230        endif ! of if (planet_type=="earth")
     231
     232        ! add random perturbation to temperature
     233        idum  = -1
     234        zz = ran1(idum)
     235        idum  = 0
     236        do l=1,llm
     237           do ij=iip2,ip1jm
     238              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
     239           enddo
     240        enddo
     241
     242        ! maintain periodicity in longitude
     243        do l=1,llm
     244           do ij=1,ip1jmp1,iip1
     245              teta(ij+iim,l)=teta(ij,l)
     246           enddo
     247        enddo
     248
     249     ENDIF ! of IF (.NOT. read_start)
     250  endif academic_case
     251
     252END SUBROUTINE iniacademic

  • LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/disvert.F90

    r1463 r1472  
    1 !
    21! $Id$
    3 !
    4       SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
    52
    6 c    Auteur :  P. Le Van .
    7 c
    8       IMPLICIT NONE
     3SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
    94
    10 #include "dimensions.h"
    11 #include "paramet.h"
    12 #include "iniprint.h"
    13 #include "logic.h"
    14 c
    15 c=======================================================================
    16 c
    17 c
    18 c    s = sigma ** kappa   :  coordonnee  verticale
    19 c    dsig(l)            : epaisseur de la couche l ds la coord.  s
    20 c    sig(l)             : sigma a l'interface des couches l et l-1
    21 c    ds(l)              : distance entre les couches l et l-1 en coord.s
    22 c
    23 c=======================================================================
    24 c
    25       REAL pa,preff
    26       REAL ap(llmp1),bp(llmp1),dpres(llm),nivsigs(llm),nivsig(llmp1)
    27       REAL presnivs(llm)
    28 c
    29 c   declarations:
    30 c   -------------
    31 c
    32       REAL sig(llm+1),dsig(llm)
    33        real zzz(1:llm+1)
    34        real dzz(1:llm)
    35       real zk,zkm1,dzk1,dzk2,k0,k1
    36 c
    37       INTEGER l
    38       REAL snorm,dsigmin
    39       REAL alpha,beta,gama,delta,deltaz,h
    40       INTEGER np,ierr
    41       REAL pi,x
     5  ! Auteur : P. Le Van
    426
    43       REAL SSUM
    44 c
    45 c-----------------------------------------------------------------------
    46 c
    47       pi=2.*ASIN(1.)
     7  IMPLICIT NONE
    488
    49       OPEN(99,file='sigma.def',status='old',form='formatted',
    50      s   iostat=ierr)
     9  include "dimensions.h"
     10  include "paramet.h"
     11  include "iniprint.h"
     12  include "logic.h"
    5113
    52 c-----------------------------------------------------------------------
    53 c   cas 1 on lit les options dans sigma.def:
    54 c   ----------------------------------------
     14  ! s = sigma ** kappa : coordonnee verticale
     15  ! dsig(l) : epaisseur de la couche l ds la coord. s
     16  ! sig(l) : sigma a l'interface des couches l et l-1
     17  ! ds(l) : distance entre les couches l et l-1 en coord.s
    5518
    56       IF (ierr.eq.0) THEN
     19  REAL pa, preff
     20  REAL ap(llmp1), bp(llmp1), dpres(llm), nivsigs(llm), nivsig(llmp1)
     21  REAL presnivs(llm)
    5722
    58       READ(99,*) h           ! hauteur d'echelle 8.
    59       READ(99,*) deltaz      ! epaiseur de la premiere couche 0.04
    60       READ(99,*) beta        ! facteur d'acroissement en haut 1.3
    61       READ(99,*) k0          ! nombre de couches dans la transition surf
    62       READ(99,*) k1          ! nombre de couches dans la transition haute
    63       CLOSE(99)
    64       alpha=deltaz/(llm*h)
    65       write(lunout,*)'h,alpha,k0,k1,beta'
     23  REAL sig(llm+1), dsig(llm)
     24  real zk, zkm1, dzk1, dzk2, k0, k1
    6625
    67 c     read(*,*) h,deltaz,beta,k0,k1 ! 8 0.04 4 20 1.2
     26  INTEGER l
     27  REAL dsigmin
     28  REAL alpha, beta, deltaz, h
     29  INTEGER iostat
     30  REAL pi, x
    6831
    69       alpha=deltaz/tanh(1./k0)*2.
    70       zkm1=0.
    71       sig(1)=1.
    72       do l=1,llm
    73         sig(l+1)=(cosh(l/k0))**(-alpha*k0/h)
    74      + *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
     32  !-----------------------------------------------------------------------
     33
     34  pi = 2 * ASIN(1.)
     35
     36  OPEN(99, file='sigma.def', status='old', form='formatted', iostat=iostat)
     37
     38  IF (iostat == 0) THEN
     39     ! cas 1 on lit les options dans sigma.def:
     40     READ(99, *) h ! hauteur d'echelle 8.
     41     READ(99, *) deltaz ! epaiseur de la premiere couche 0.04
     42     READ(99, *) beta ! facteur d'acroissement en haut 1.3
     43     READ(99, *) k0 ! nombre de couches dans la transition surf
     44     READ(99, *) k1 ! nombre de couches dans la transition haute
     45     CLOSE(99)
     46     alpha=deltaz/(llm*h)
     47     write(lunout, *)'h, alpha, k0, k1, beta'
     48
     49     alpha=deltaz/tanh(1./k0)*2.
     50     zkm1=0.
     51     sig(1)=1.
     52     do l=1, llm
     53        sig(l+1)=(cosh(l/k0))**(-alpha*k0/h) &
     54             *exp(-alpha/h*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta))
    7555        zk=-h*log(sig(l+1))
    7656
    7757        dzk1=alpha*tanh(l/k0)
    7858        dzk2=alpha*tanh((llm-k1)/k0)*beta**(l-(llm-k1))/log(beta)
    79         write(lunout,*)l,sig(l+1),zk,zk-zkm1,dzk1,dzk2
     59        write(lunout, *)l, sig(l+1), zk, zk-zkm1, dzk1, dzk2
    8060        zkm1=zk
    81       enddo
     61     enddo
    8262
    83       sig(llm+1)=0.
     63     sig(llm+1)=0.
    8464
    85 c
    86        DO 2  l = 1, llm
    87        dsig(l) = sig(l)-sig(l+1)
    88    2   CONTINUE
    89 c
     65     DO l = 1, llm
     66        dsig(l) = sig(l)-sig(l+1)
     67     end DO
     68  ELSE
     69     if (ok_strato) then
     70        if (llm==39) then
     71           dsigmin=0.3
     72        else if (llm==50) then
     73           dsigmin=1.
     74        else
     75           WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
     76           dsigmin=1.
     77        endif
     78        WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
     79     endif
    9080
    91       ELSE
    92 c-----------------------------------------------------------------------
    93 c   cas 2 ancienne discretisation (LMD5...):
    94 c   ----------------------------------------
     81     DO l = 1, llm
     82        x = pi * (l - 0.5) / (llm + 1)
    9583
    96       WRITE(LUNOUT,*)'WARNING!!! Ancienne discretisation verticale'
     84        IF (ok_strato) THEN
     85           dsig(l) =(dsigmin + 7. * SIN(x)**2) &
     86                *(1. - tanh(2 * x / pi - 1.))**2 / 4.
     87        ELSE
     88           dsig(l) = 1.0 + 7.0 * SIN(x)**2
     89        ENDIF
     90     ENDDO
     91     dsig = dsig / sum(dsig)
     92     sig(llm+1) = 0.
     93     DO l = llm, 1, -1
     94        sig(l) = sig(l+1) + dsig(l)
     95     ENDDO
     96  ENDIF
    9797
    98       if (ok_strato) then
    99          if (llm==39) then
    100             dsigmin=0.3
    101          else if (llm==50) then
    102             dsigmin=1.
    103          else
    104             WRITE(LUNOUT,*) 'ATTENTION discretisation z a ajuster'
    105             dsigmin=1.
    106          endif
    107          WRITE(LUNOUT,*) 'Discretisation verticale DSIGMIN=',dsigmin
    108       endif
     98  DO l=1, llm
     99     nivsigs(l) = REAL(l)
     100  ENDDO
    109101
    110       h=7.
    111       snorm  = 0.
    112       DO l = 1, llm
    113          x = 2.*asin(1.) * (REAL(l)-0.5) / REAL(llm+1)
     102  DO l=1, llmp1
     103     nivsig(l)= REAL(l)
     104  ENDDO
    114105
    115          IF (ok_strato) THEN
    116            dsig(l) =(dsigmin + 7.0 * SIN(x)**2)
    117      &            *(0.5*(1.-tanh(1.*(x-asin(1.))/asin(1.))))**2       
    118          ELSE
    119            dsig(l) = 1.0 + 7.0 * SIN(x)**2
    120          ENDIF
     106  ! .... Calculs de ap(l) et de bp(l) ....
     107  ! ..... pa et preff sont lus sur les fichiers start par lectba .....
    121108
    122          snorm = snorm + dsig(l)
    123       ENDDO
    124       snorm = 1./snorm
    125       DO l = 1, llm
    126          dsig(l) = dsig(l)*snorm
    127       ENDDO
    128       sig(llm+1) = 0.
    129       DO l = llm, 1, -1
    130          sig(l) = sig(l+1) + dsig(l)
    131       ENDDO
     109  bp(llmp1) = 0.
    132110
    133       ENDIF
     111  DO l = 1, llm
     112     bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
     113     ap(l) = pa * ( sig(l) - bp(l) )
     114  ENDDO
    134115
     116  bp(1)=1.
     117  ap(1)=0.
    135118
    136       DO l=1,llm
    137         nivsigs(l) = REAL(l)
    138       ENDDO
     119  ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
    139120
    140       DO l=1,llmp1
    141         nivsig(l)= REAL(l)
    142       ENDDO
     121  write(lunout, *)' BP '
     122  write(lunout, *) bp
     123  write(lunout, *)' AP '
     124  write(lunout, *) ap
    143125
    144 c
    145 c    ....  Calculs  de ap(l) et de bp(l)  ....
    146 c    .........................................
    147 c
    148 c
    149 c   .....  pa et preff sont lus  sur les fichiers start par lectba  .....
    150 c
     126  write(lunout, *) 'Niveaux de pressions approximatifs aux centres des'
     127  write(lunout, *)'couches calcules pour une pression de surface =', preff
     128  write(lunout, *) 'et altitudes equivalentes pour une hauteur d echelle de'
     129  write(lunout, *)'8km'
     130  DO l = 1, llm
     131     dpres(l) = bp(l) - bp(l+1)
     132     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
     133     write(lunout, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
     134          log(preff/presnivs(l))*8. &
     135          , ' DZ ~ ', 8.*log((ap(l)+bp(l)*preff)/ &
     136          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
     137  ENDDO
    151138
    152       bp(llmp1) =   0.
     139  write(lunout, *) 'PRESNIVS '
     140  write(lunout, *) presnivs
    153141
    154       DO l = 1, llm
    155 cc
    156 ccc    ap(l) = 0.
    157 ccc    bp(l) = sig(l)
    158 
    159       bp(l) = EXP( 1. -1./( sig(l)*sig(l)) )
    160       ap(l) = pa * ( sig(l) - bp(l) )
    161 c
    162       ENDDO
    163 
    164       bp(1)=1.
    165       ap(1)=0.
    166 
    167       ap(llmp1) = pa * ( sig(llmp1) - bp(llmp1) )
    168 
    169       write(lunout,*)' BP '
    170       write(lunout,*)  bp
    171       write(lunout,*)' AP '
    172       write(lunout,*)  ap
    173 
    174       write(lunout,*)
    175      .'Niveaux de pressions approximatifs aux centres des'
    176       write(lunout,*)'couches calcules pour une pression de surface =',
    177      .                 preff
    178       write(lunout,*)
    179      .     'et altitudes equivalentes pour une hauteur d echelle de'
    180       write(lunout,*)'8km'
    181       DO l = 1, llm
    182        dpres(l) = bp(l) - bp(l+1)
    183        presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
    184        write(lunout,*)'PRESNIVS(',l,')=',presnivs(l),'    Z ~ ',
    185      .        log(preff/presnivs(l))*8.
    186      .  ,'   DZ ~ ',8.*log((ap(l)+bp(l)*preff)/
    187      .       max(ap(l+1)+bp(l+1)*preff,1.e-10))
    188       ENDDO
    189 
    190       write(lunout,*)' PRESNIVS '
    191       write(lunout,*)presnivs
    192 
    193       RETURN
    194       END
     142END SUBROUTINE disvert

  • LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3dpar/iniacademic.F90

    r1463 r1472  
    22! $Id$
    33!
    4 c
    5 c
    6       SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
    7 
    8       USE filtreg_mod
    9       USE infotrac, ONLY : nqtot
    10       USE control_mod, ONLY: day_step,planet_type
     4SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
     5
     6  USE filtreg_mod
     7  USE infotrac, ONLY : nqtot
     8  USE control_mod, ONLY: day_step,planet_type
    119#ifdef CPP_IOIPSL
    12       USE IOIPSL
     10  USE IOIPSL
    1311#else
    14 ! if not using IOIPSL, we still need to use (a local version of) getin
    15       USE ioipsl_getincom
     12  ! if not using IOIPSL, we still need to use (a local version of) getin
     13  USE ioipsl_getincom
    1614#endif
    17       USE Write_Field
    18  
    19 
    20 c%W%    %G%
    21 c=======================================================================
    22 c
    23 c   Author:    Frederic Hourdin      original: 15/01/93
    24 c   -------
    25 c
    26 c   Subject:
    27 c   ------
    28 c
    29 c   Method:
    30 c   --------
    31 c
    32 c   Interface:
    33 c   ----------
    34 c
    35 c      Input:
    36 c      ------
    37 c
    38 c      Output:
    39 c      -------
    40 c
    41 c=======================================================================
    42       IMPLICIT NONE
    43 c-----------------------------------------------------------------------
    44 c   Declararations:
    45 c   ---------------
    46 
    47 #include "dimensions.h"
    48 #include "paramet.h"
    49 #include "comvert.h"
    50 #include "comconst.h"
    51 #include "comgeom.h"
    52 #include "academic.h"
    53 #include "ener.h"
    54 #include "temps.h"
    55 #include "iniprint.h"
    56 #include "logic.h"
    57 
    58 c   Arguments:
    59 c   ----------
    60 
    61       real time_0
    62 
    63 c   variables dynamiques
    64       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
    65       REAL teta(ip1jmp1,llm)                 ! temperature potentielle
    66       REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
    67       REAL ps(ip1jmp1)                       ! pression  au sol
    68       REAL masse(ip1jmp1,llm)                ! masse d'air
    69       REAL phis(ip1jmp1)                     ! geopotentiel au sol
    70 
    71 c   Local:
    72 c   ------
    73 
    74       REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
    75       REAL pks(ip1jmp1)                      ! exner au  sol
    76       REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
    77       REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
    78       REAL phi(ip1jmp1,llm)                  ! geopotentiel
    79       REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
    80       real tetajl(jjp1,llm)
    81       INTEGER i,j,l,lsup,ij
    82 
    83       REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
    84       REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
    85       LOGICAL ok_geost             ! Initialisation vent geost. ou nul
    86       LOGICAL ok_pv                ! Polar Vortex
    87       REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex
    88      
    89       real zz,ran1
    90       integer idum
    91 
    92       REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
    93 
    94 c-----------------------------------------------------------------------
    95 ! 1. Initializations for Earth-like case
    96 ! --------------------------------------
    97 c
    98         ! initialize planet radius, rotation rate,...
    99         call conf_planete
    100 
    101         time_0=0.
    102         day_ref=1
    103         annee_ref=0
    104 
    105         im         = iim
    106         jm         = jjm
    107         day_ini    = 1
    108         dtvr    = daysec/REAL(day_step)
    109         zdtvr=dtvr
    110         etot0      = 0.
    111         ptot0      = 0.
    112         ztot0      = 0.
    113         stot0      = 0.
    114         ang0       = 0.
    115 
    116         if (llm.eq.1) then
    117           ! specific initializations for the shallow water case
    118           kappa=1
     15  USE Write_Field
     16
     17  !   Author:    Frederic Hourdin      original: 15/01/93
     18
     19  IMPLICIT NONE
     20
     21  !   Declararations:
     22  !   ---------------
     23
     24  include "dimensions.h"
     25  include "paramet.h"
     26  include "comvert.h"
     27  include "comconst.h"
     28  include "comgeom.h"
     29  include "academic.h"
     30  include "ener.h"
     31  include "temps.h"
     32  include "iniprint.h"
     33  include "logic.h"
     34
     35  !   Arguments:
     36  !   ----------
     37
     38  real time_0
     39
     40  !   variables dynamiques
     41  REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
     42  REAL teta(ip1jmp1,llm)                 ! temperature potentielle
     43  REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
     44  REAL ps(ip1jmp1)                       ! pression  au sol
     45  REAL masse(ip1jmp1,llm)                ! masse d'air
     46  REAL phis(ip1jmp1)                     ! geopotentiel au sol
     47
     48  !   Local:
     49  !   ------
     50
     51  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
     52  REAL pks(ip1jmp1)                      ! exner au  sol
     53  REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
     54  REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
     55  REAL phi(ip1jmp1,llm)                  ! geopotentiel
     56  REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
     57  real tetajl(jjp1,llm)
     58  INTEGER i,j,l,lsup,ij
     59
     60  REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
     61  REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
     62  LOGICAL ok_geost             ! Initialisation vent geost. ou nul
     63  LOGICAL ok_pv                ! Polar Vortex
     64  REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex
     65
     66  real zz,ran1
     67  integer idum
     68
     69  REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
     70
     71  !-----------------------------------------------------------------------
     72  ! 1. Initializations for Earth-like case
     73  ! --------------------------------------
     74  !
     75  ! initialize planet radius, rotation rate,...
     76  call conf_planete
     77
     78  time_0=0.
     79  day_ref=1
     80  annee_ref=0
     81
     82  im         = iim
     83  jm         = jjm
     84  day_ini    = 1
     85  dtvr    = daysec/REAL(day_step)
     86  zdtvr=dtvr
     87  etot0      = 0.
     88  ptot0      = 0.
     89  ztot0      = 0.
     90  stot0      = 0.
     91  ang0       = 0.
     92
     93  if (llm == 1) then
     94     ! specific initializations for the shallow water case
     95     kappa=1
     96  endif
     97
     98  CALL iniconst
     99  CALL inigeom
     100  CALL inifilr
     101
     102  if (llm == 1) then
     103     ! initialize fields for the shallow water case, if required
     104     if (.not.read_start) then
     105        phis(:)=0.
     106        q(:,:,:)=0
     107        CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
     108     endif
     109  endif
     110
     111  academic_case: if (iflag_phys == 2) then
     112     ! initializations
     113
     114     ! 1. local parameters
     115     ! by convention, winter is in the southern hemisphere
     116     ! Geostrophic wind or no wind?
     117     ok_geost=.TRUE.
     118     CALL getin('ok_geost',ok_geost)
     119     ! Constants for Newtonian relaxation and friction
     120     k_f=1.                !friction
     121     CALL getin('k_j',k_f)
     122     k_f=1./(daysec*k_f)
     123     k_c_s=4.  !cooling surface
     124     CALL getin('k_c_s',k_c_s)
     125     k_c_s=1./(daysec*k_c_s)
     126     k_c_a=40. !cooling free atm
     127     CALL getin('k_c_a',k_c_a)
     128     k_c_a=1./(daysec*k_c_a)
     129     ! Constants for Teta equilibrium profile
     130     teta0=315.     ! mean Teta (S.H. 315K)
     131     CALL getin('teta0',teta0)
     132     ttp=200.       ! Tropopause temperature (S.H. 200K)
     133     CALL getin('ttp',ttp)
     134     eps=0.         ! Deviation to N-S symmetry(~0-20K)
     135     CALL getin('eps',eps)
     136     delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
     137     CALL getin('delt_y',delt_y)
     138     delt_z=10.     ! Vertical Gradient (S.H. 10K)
     139     CALL getin('delt_z',delt_z)
     140     ! Polar vortex
     141     ok_pv=.false.
     142     CALL getin('ok_pv',ok_pv)
     143     phi_pv=-50.            ! Latitude of edge of vortex
     144     CALL getin('phi_pv',phi_pv)
     145     phi_pv=phi_pv*pi/180.
     146     dphi_pv=5.             ! Width of the edge
     147     CALL getin('dphi_pv',dphi_pv)
     148     dphi_pv=dphi_pv*pi/180.
     149     gam_pv=4.              ! -dT/dz vortex (in K/km)
     150     CALL getin('gam_pv',gam_pv)
     151
     152     ! 2. Initialize fields towards which to relax
     153     ! Friction
     154     knewt_g=k_c_a
     155     DO l=1,llm
     156        zsig=presnivs(l)/preff
     157        knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
     158        kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
     159     ENDDO
     160     DO j=1,jjp1
     161        clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
     162     ENDDO
     163
     164     ! Potential temperature
     165     DO l=1,llm
     166        zsig=presnivs(l)/preff
     167        tetastrat=ttp*zsig**(-kappa)
     168        tetapv=tetastrat
     169        IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
     170           tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
     171        ENDIF
     172        DO j=1,jjp1
     173           ! Troposphere
     174           ddsin=sin(rlatu(j))
     175           tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
     176                -delt_z*(1.-ddsin*ddsin)*log(zsig)
     177           ! Profil stratospherique isotherme (+vortex)
     178           w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
     179           tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
     180           tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
     181        ENDDO
     182     ENDDO
     183
     184     !          CALL writefield('theta_eq',tetajl)
     185
     186     do l=1,llm
     187        do j=1,jjp1
     188           do i=1,iip1
     189              ij=(j-1)*iip1+i
     190              tetarappel(ij,l)=tetajl(j,l)
     191           enddo
     192        enddo
     193     enddo
     194
     195     ! 3. Initialize fields (if necessary)
     196     IF (.NOT. read_start) THEN
     197        ! surface pressure
     198        ps(:)=preff
     199        ! ground geopotential
     200        phis(:)=0.
     201
     202        CALL pression ( ip1jmp1, ap, bp, ps, p       )
     203        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     204        CALL massdair(p,masse)
     205
     206        ! bulk initialization of temperature
     207        teta(:,:)=tetarappel(:,:)
     208
     209        ! geopotential
     210        CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
     211
     212        ! winds
     213        if (ok_geost) then
     214           call ugeostr(phi,ucov)
     215        else
     216           ucov(:,:)=0.
    119217        endif
    120        
    121         CALL iniconst
    122         CALL inigeom
    123         CALL inifilr
    124 
    125         if (llm.eq.1) then
    126           ! initialize fields for the shallow water case, if required
    127           if (.not.read_start) then
    128             phis(:)=0.
    129             q(:,:,:)=0
    130             CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
    131           endif
    132         endif
    133 
    134         if (iflag_phys.eq.2) then
    135           ! initializations for the academic case
    136          
    137 !         if (planet_type=="earth") then
    138 
    139           ! 1. local parameters
    140           ! by convention, winter is in the southern hemisphere
    141           ! Geostrophic wind or no wind?
    142           ok_geost=.TRUE.
    143           CALL getin('ok_geost',ok_geost)
    144           ! Constants for Newtonian relaxation and friction
    145           k_f=1.                !friction
    146           CALL getin('k_j',k_f)
    147           k_f=1./(daysec*k_f)
    148           k_c_s=4.  !cooling surface
    149           CALL getin('k_c_s',k_c_s)
    150           k_c_s=1./(daysec*k_c_s)
    151           k_c_a=40. !cooling free atm
    152           CALL getin('k_c_a',k_c_a)
    153           k_c_a=1./(daysec*k_c_a)
    154           ! Constants for Teta equilibrium profile
    155           teta0=315.     ! mean Teta (S.H. 315K)
    156           CALL getin('teta0',teta0)
    157           ttp=200.       ! Tropopause temperature (S.H. 200K)
    158           CALL getin('ttp',ttp)
    159           eps=0.         ! Deviation to N-S symmetry(~0-20K)
    160           CALL getin('eps',eps)
    161           delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
    162           CALL getin('delt_y',delt_y)
    163           delt_z=10.     ! Vertical Gradient (S.H. 10K)
    164           CALL getin('delt_z',delt_z)
    165           ! Polar vortex
    166           ok_pv=.false.
    167           CALL getin('ok_pv',ok_pv)
    168           phi_pv=-50.            ! Latitude of edge of vortex
    169           CALL getin('phi_pv',phi_pv)
    170           phi_pv=phi_pv*pi/180.
    171           dphi_pv=5.             ! Width of the edge
    172           CALL getin('dphi_pv',dphi_pv)
    173           dphi_pv=dphi_pv*pi/180.
    174           gam_pv=4.              ! -dT/dz vortex (in K/km)
    175           CALL getin('gam_pv',gam_pv)
    176          
    177           ! 2. Initialize fields towards which to relax
    178           ! Friction
    179           knewt_g=k_c_a
    180           DO l=1,llm
    181            zsig=presnivs(l)/preff
    182            knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
    183            kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
    184           ENDDO
    185           DO j=1,jjp1
    186             clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
    187           ENDDO
    188          
    189           ! Potential temperature
    190           DO l=1,llm
    191            zsig=presnivs(l)/preff
    192            tetastrat=ttp*zsig**(-kappa)
    193            tetapv=tetastrat
    194            IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
    195                tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
    196            ENDIF
    197            DO j=1,jjp1
    198              ! Troposphere
    199              ddsin=sin(rlatu(j))
    200              tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin
    201      &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
    202              ! Profil stratospherique isotherme (+vortex)
    203              w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
    204              tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
    205              tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
    206            ENDDO
    207           ENDDO ! of DO l=1,llm
    208  
    209 !          CALL writefield('theta_eq',tetajl)
    210 
    211           do l=1,llm
    212             do j=1,jjp1
    213               do i=1,iip1
    214                  ij=(j-1)*iip1+i
    215                  tetarappel(ij,l)=tetajl(j,l)
    216               enddo
    217             enddo
    218           enddo
    219 
    220 
    221 !         else
    222 !          write(lunout,*)"iniacademic: planet types other than earth",
    223 !     &                   " not implemented (yet)."
    224 !          stop
    225 !         endif ! of if (planet_type=="earth")
    226 
    227           ! 3. Initialize fields (if necessary)
    228           IF (.NOT. read_start) THEN
    229             ! surface pressure
    230             ps(:)=preff
    231             ! ground geopotential
    232             phis(:)=0.
    233            
    234             CALL pression ( ip1jmp1, ap, bp, ps, p       )
    235             CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    236             CALL massdair(p,masse)
    237 
    238             ! bulk initialization of temperature
    239             teta(:,:)=tetarappel(:,:)
    240            
    241             ! geopotential
    242             CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
    243            
    244             ! winds
    245             if (ok_geost) then
    246               call ugeostr(phi,ucov)
    247             else
    248               ucov(:,:)=0.
    249             endif
    250             vcov(:,:)=0.
    251            
    252             ! bulk initialization of tracers
    253             if (planet_type=="earth") then
    254               ! Earth: first two tracers will be water
    255               do i=1,nqtot
    256                 if (i.eq.1) q(:,:,i)=1.e-10
    257                 if (i.eq.2) q(:,:,i)=1.e-15
    258                 if (i.gt.2) q(:,:,i)=0.
    259               enddo
    260             else
    261               q(:,:,:)=0
    262             endif ! of if (planet_type=="earth")
    263 
    264             ! add random perturbation to temperature
    265             idum  = -1
    266             zz = ran1(idum)
    267             idum  = 0
    268             do l=1,llm
    269               do ij=iip2,ip1jm
    270                 teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
    271               enddo
    272             enddo
    273 
    274             ! maintain periodicity in longitude
    275             do l=1,llm
    276               do ij=1,ip1jmp1,iip1
    277                 teta(ij+iim,l)=teta(ij,l)
    278               enddo
    279             enddo
    280 
    281           ENDIF ! of IF (.NOT. read_start)
    282         endif ! of if (iflag_phys.eq.2)
    283        
    284       END
    285 c-----------------------------------------------------------------------
     218        vcov(:,:)=0.
     219
     220        ! bulk initialization of tracers
     221        if (planet_type=="earth") then
     222           ! Earth: first two tracers will be water
     223           do i=1,nqtot
     224              if (i == 1) q(:,:,i)=1.e-10
     225              if (i == 2) q(:,:,i)=1.e-15
     226              if (i.gt.2) q(:,:,i)=0.
     227           enddo
     228        else
     229           q(:,:,:)=0
     230        endif ! of if (planet_type=="earth")
     231
     232        ! add random perturbation to temperature
     233        idum  = -1
     234        zz = ran1(idum)
     235        idum  = 0
     236        do l=1,llm
     237           do ij=iip2,ip1jm
     238              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
     239           enddo
     240        enddo
     241
     242        ! maintain periodicity in longitude
     243        do l=1,llm
     244           do ij=1,ip1jmp1,iip1
     245              teta(ij+iim,l)=teta(ij,l)
     246           enddo
     247        enddo
     248
     249     ENDIF ! of IF (.NOT. read_start)
     250  endif academic_case
     251
     252END SUBROUTINE iniacademic

  • LMDZ5/branches/LMDZ5V2.0-dev/libf/phylmd/fisrtilp.F90

    r1463 r1472  
    22! $Id$
    33!
    4 c
    5       SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs,
    6      s                   d_t, d_q, d_ql, rneb, radliq, rain, snow,
    7      s                   pfrac_impa, pfrac_nucl, pfrac_1nucl,
    8      s                   frac_impa, frac_nucl,
    9      s                   prfl, psfl, rhcl, zqta, fraca,
    10      s                   ztv, zpspsk, ztla, zthl, iflag_cldcon)
    11 
    12 c
    13       USE dimphy
    14       IMPLICIT none
    15 c======================================================================
    16 c Auteur(s): Z.X. Li (LMD/CNRS)
    17 c Date: le 20 mars 1995
    18 c Objet: condensation et precipitation stratiforme.
    19 c        schema de nuage
    20 c======================================================================
    21 c======================================================================
    22 cym#include "dimensions.h"
    23 cym#include "dimphy.h"
    24 #include "YOMCST.h"
    25 #include "tracstoke.h"
    26 #include "fisrtilp.h"
    27 c
    28 c Arguments:
    29 c
    30       REAL dtime ! intervalle du temps (s)
    31       REAL paprs(klon,klev+1) ! pression a inter-couche
    32       REAL pplay(klon,klev) ! pression au milieu de couche
    33       REAL t(klon,klev) ! temperature (K)
    34       REAL q(klon,klev) ! humidite specifique (kg/kg)
    35       REAL d_t(klon,klev) ! incrementation de la temperature (K)
    36       REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
    37       REAL d_ql(klon,klev) ! incrementation de l'eau liquide
    38       REAL rneb(klon,klev) ! fraction nuageuse
    39       REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
    40       REAL rhcl(klon,klev) ! humidite relative en ciel clair
    41       REAL rain(klon) ! pluies (mm/s)
    42       REAL snow(klon) ! neige (mm/s)
    43       REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
    44       REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
    45       REAL ztv(klon,klev)
    46       REAL zqta(klon,klev),fraca(klon,klev)
    47       REAL sigma1(klon,klev),sigma2(klon,klev)
    48       REAL qltot(klon,klev),ctot(klon,klev)
    49       REAL zpspsk(klon,klev),ztla(klon,klev)
    50       REAL zthl(klon,klev)
    51 
    52       logical lognormale(klon)
    53 
    54 cAA
    55 c Coeffients de fraction lessivee : pour OFF-LINE
    56 c
    57       REAL pfrac_nucl(klon,klev)
    58       REAL pfrac_1nucl(klon,klev)
    59       REAL pfrac_impa(klon,klev)
    60 c
    61 c Fraction d'aerosols lessivee par impaction et par nucleation
    62 c POur ON-LINE
    63 c
    64       REAL frac_impa(klon,klev)
    65       REAL frac_nucl(klon,klev)
    66       real zct      ,zcl
    67 cAA
    68 c
    69 c Options du programme:
    70 c
    71       REAL seuil_neb ! un nuage existe vraiment au-dela
    72       PARAMETER (seuil_neb=0.001)
    73 
    74       INTEGER ninter ! sous-intervals pour la precipitation
    75       INTEGER ncoreczq 
    76       INTEGER iflag_cldcon
    77       PARAMETER (ninter=5)
    78       LOGICAL evap_prec ! evaporation de la pluie
    79       PARAMETER (evap_prec=.TRUE.)
    80       REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
    81       logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
    82 
    83       real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
    84       real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
    85       real erf   
    86       REAL qcloud(klon)
    87 c
    88       LOGICAL cpartiel ! condensation partielle
    89       PARAMETER (cpartiel=.TRUE.)
    90       REAL t_coup
    91       PARAMETER (t_coup=234.0)
    92 c
    93 c Variables locales:
    94 c
    95       INTEGER i, k, n, kk
    96       REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5   
    97       REAL zrfl(klon), zrfln(klon), zqev, zqevt
    98       REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
    99       REAL ztglace, zt(klon)
    100       INTEGER nexpo ! exponentiel pour glace/eau
    101       REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
    102       REAL zchau      ,zfroi      ,zfice(klon),zneb(klon)
    103 c
    104       LOGICAL appel1er
    105       SAVE appel1er
    106 c$OMP THREADPRIVATE(appel1er)
    107 c
    108 c---------------------------------------------------------------
    109 c
    110 cAA Variables traceurs:
    111 cAA  Provisoire !!! Parametres alpha du lessivage
    112 cAA  A priori on a 4 scavenging # possibles
    113 c
    114       REAL a_tr_sca(4)
    115       save a_tr_sca
    116 c$OMP THREADPRIVATE(a_tr_sca)
    117 c
    118 c Variables intermediaires
    119 c
    120       REAL zalpha_tr
    121       REAL zfrac_lessi
    122       REAL zprec_cond(klon)
    123 cAA
    124       REAL zmair, zcpair, zcpeau
    125 C     Pour la conversion eau-neige
    126       REAL zlh_solid(klon), zm_solid
    127 cIM
    128 cym      INTEGER klevm1
    129 c---------------------------------------------------------------
    130 c
    131 c Fonctions en ligne:
    132 c
    133       REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
    134       REAL zzz
    135 #include "YOETHF.h"
    136 #include "FCTTRE.h"
    137       fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
    138       fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
    139 c
    140       DATA appel1er /.TRUE./
    141 cym
    142       zdelq=0.0
    143      
    144       print*,'NUAGES4 A. JAM'
    145       IF (appel1er) THEN
    146 c
    147          PRINT*, 'fisrtilp, ninter:', ninter
    148          PRINT*, 'fisrtilp, evap_prec:', evap_prec
    149          PRINT*, 'fisrtilp, cpartiel:', cpartiel
    150          IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
    151           PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
    152           PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
    153 c         CALL abort
    154          ENDIF
    155          appel1er = .FALSE.
    156 c
    157 cAA initialiation provisoire
    158        a_tr_sca(1) = -0.5
    159        a_tr_sca(2) = -0.5
    160        a_tr_sca(3) = -0.5
    161        a_tr_sca(4) = -0.5
    162 c
    163 cAA Initialisation a 1 des coefs des fractions lessivees
    164 c
    165 !cdir collapse
    166       DO k = 1, klev
    167        DO i = 1, klon
    168           pfrac_nucl(i,k)=1.
    169           pfrac_1nucl(i,k)=1.
    170           pfrac_impa(i,k)=1.
    171        ENDDO
    172       ENDDO
    173 
    174       ENDIF          !  test sur appel1er
    175 c
    176 cMAf Initialisation a 0 de zoliq
    177 c      DO i = 1, klon
    178 c         zoliq(i)=0.
    179 c      ENDDO
    180 c Determiner les nuages froids par leur temperature
    181 c  nexpo regle la raideur de la transition eau liquide / eau glace.
    182 c
    183       ztglace = RTT - 15.0
    184       nexpo = 6
    185 ccc      nexpo = 1
    186 c
    187 c Initialiser les sorties:
    188 c
    189 !cdir collapse
    190       DO k = 1, klev+1
    191       DO i = 1, klon
    192          prfl(i,k) = 0.0
    193          psfl(i,k) = 0.0
    194       ENDDO
    195       ENDDO
    196 
    197 !cdir collapse
    198       DO k = 1, klev
    199       DO i = 1, klon
    200          d_t(i,k) = 0.0
    201          d_q(i,k) = 0.0
    202          d_ql(i,k) = 0.0
    203          rneb(i,k) = 0.0
    204          radliq(i,k) = 0.0
    205          frac_nucl(i,k) = 1.
    206          frac_impa(i,k) = 1.
    207       ENDDO
    208       ENDDO
    209       DO i = 1, klon
    210          rain(i) = 0.0
    211          snow(i) = 0.0
    212          zoliq(i)=0.
    213 c     ENDDO
    214 c
    215 c Initialiser le flux de precipitation a zero
    216 c
    217 c     DO i = 1, klon
    218          zrfl(i) = 0.0
    219          zneb(i) = seuil_neb
    220       ENDDO
    221 c
    222 c
    223 cAA Pour plus de securite
    224 
    225       zalpha_tr   = 0.
    226       zfrac_lessi = 0.
    227 
    228 cAA----------------------------------------------------------
    229 c
    230       ncoreczq=0
    231 c Boucle verticale (du haut vers le bas)
    232 c
    233 cIM : klevm1
    234 cym      klevm1=klev-1
    235       DO 9999 k = klev, 1, -1
    236 c
    237 cAA----------------------------------------------------------
    238 c
    239       DO i = 1, klon
    240          zt(i)=t(i,k)
    241          zq(i)=q(i,k)
    242       ENDDO
    243 c
    244 c Calculer la varition de temp. de l'air du a la chaleur sensible
    245 C transporter par la pluie.
    246 C Il resterait a rajouter cet effet de la chaleur sensible sur les
    247 C flux de surface, du a la diff. de temp. entre le 1er niveau et la
    248 C surface.
    249 C
    250       IF(k.LE.klevm1) THEN         
    251          DO i = 1, klon
    252 cIM
    253             zmair=(paprs(i,k)-paprs(i,k+1))/RG
    254             zcpair=RCPD*(1.0+RVTMP2*zq(i))
    255             zcpeau=RCPD*RVTMP2
    256             zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau
    257      $           + zmair*zcpair*zt(i) )
    258      $           / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
    259 C     C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
    260          ENDDO
    261       ENDIF
    262 c
    263 c
    264 c Calculer l'evaporation de la precipitation
    265 c
    266 
    267 
    268       IF (evap_prec) THEN
    269       DO i = 1, klon
    270       IF (zrfl(i) .GT.0.) THEN
    271          IF (thermcep) THEN
    272            zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
    273            zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
    274            zqs(i)=MIN(0.5,zqs(i))
    275            zcor=1./(1.-RETV*zqs(i))
    276            zqs(i)=zqs(i)*zcor
    277          ELSE
    278            IF (zt(i) .LT. t_coup) THEN
    279               zqs(i) = qsats(zt(i)) / pplay(i,k)
    280            ELSE
    281               zqs(i) = qsatl(zt(i)) / pplay(i,k)
     4!
     5SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
     6     d_t, d_q, d_ql, rneb, radliq, rain, snow, &
     7     pfrac_impa, pfrac_nucl, pfrac_1nucl, &
     8     frac_impa, frac_nucl, &
     9     prfl, psfl, rhcl, zqta, fraca, &
     10     ztv, zpspsk, ztla, zthl, iflag_cldcon)
     11
     12  !
     13  USE dimphy
     14  IMPLICIT none
     15  !======================================================================
     16  ! Auteur(s): Z.X. Li (LMD/CNRS)
     17  ! Date: le 20 mars 1995
     18  ! Objet: condensation et precipitation stratiforme.
     19  !        schema de nuage
     20  !======================================================================
     21  !======================================================================
     22  !ym include "dimensions.h"
     23  !ym include "dimphy.h"
     24  include "YOMCST.h"
     25  include "tracstoke.h"
     26  include "fisrtilp.h"
     27  !
     28  ! Arguments:
     29  !
     30  REAL dtime ! intervalle du temps (s)
     31  REAL paprs(klon,klev+1) ! pression a inter-couche
     32  REAL pplay(klon,klev) ! pression au milieu de couche
     33  REAL t(klon,klev) ! temperature (K)
     34  REAL q(klon,klev) ! humidite specifique (kg/kg)
     35  REAL d_t(klon,klev) ! incrementation de la temperature (K)
     36  REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
     37  REAL d_ql(klon,klev) ! incrementation de l'eau liquide
     38  REAL rneb(klon,klev) ! fraction nuageuse
     39  REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
     40  REAL rhcl(klon,klev) ! humidite relative en ciel clair
     41  REAL rain(klon) ! pluies (mm/s)
     42  REAL snow(klon) ! neige (mm/s)
     43  REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
     44  REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
     45  REAL ztv(klon,klev)
     46  REAL zqta(klon,klev),fraca(klon,klev)
     47  REAL sigma1(klon,klev),sigma2(klon,klev)
     48  REAL qltot(klon,klev),ctot(klon,klev)
     49  REAL zpspsk(klon,klev),ztla(klon,klev)
     50  REAL zthl(klon,klev)
     51
     52  logical lognormale(klon)
     53
     54  !AA
     55  ! Coeffients de fraction lessivee : pour OFF-LINE
     56  !
     57  REAL pfrac_nucl(klon,klev)
     58  REAL pfrac_1nucl(klon,klev)
     59  REAL pfrac_impa(klon,klev)
     60  !
     61  ! Fraction d'aerosols lessivee par impaction et par nucleation
     62  ! POur ON-LINE
     63  !
     64  REAL frac_impa(klon,klev)
     65  REAL frac_nucl(klon,klev)
     66  real zct      ,zcl
     67  !AA
     68  !
     69  ! Options du programme:
     70  !
     71  REAL seuil_neb ! un nuage existe vraiment au-dela
     72  PARAMETER (seuil_neb=0.001)
     73
     74  INTEGER ninter ! sous-intervals pour la precipitation
     75  INTEGER ncoreczq 
     76  INTEGER iflag_cldcon
     77  PARAMETER (ninter=5)
     78  LOGICAL evap_prec ! evaporation de la pluie
     79  PARAMETER (evap_prec=.TRUE.)
     80  REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
     81  logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
     82
     83  real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
     84  real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
     85  real erf   
     86  REAL qcloud(klon)
     87  !
     88  LOGICAL cpartiel ! condensation partielle
     89  PARAMETER (cpartiel=.TRUE.)
     90  REAL t_coup
     91  PARAMETER (t_coup=234.0)
     92  !
     93  ! Variables locales:
     94  !
     95  INTEGER i, k, n, kk
     96  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5   
     97  REAL zrfl(klon), zrfln(klon), zqev, zqevt
     98  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
     99  REAL ztglace, zt(klon)
     100  INTEGER nexpo ! exponentiel pour glace/eau
     101  REAL zdz(klon),zrho(klon),ztot      , zrhol(klon)
     102  REAL zchau      ,zfroi      ,zfice(klon),zneb(klon)
     103  !
     104  LOGICAL appel1er
     105  SAVE appel1er
     106  !$OMP THREADPRIVATE(appel1er)
     107  !
     108  !---------------------------------------------------------------
     109  !
     110  !AA Variables traceurs:
     111  !AA  Provisoire !!! Parametres alpha du lessivage
     112  !AA  A priori on a 4 scavenging # possibles
     113  !
     114  REAL a_tr_sca(4)
     115  save a_tr_sca
     116  !$OMP THREADPRIVATE(a_tr_sca)
     117  !
     118  ! Variables intermediaires
     119  !
     120  REAL zalpha_tr
     121  REAL zfrac_lessi
     122  REAL zprec_cond(klon)
     123  !AA
     124  REAL zmair, zcpair, zcpeau
     125  !     Pour la conversion eau-neige
     126  REAL zlh_solid(klon), zm_solid
     127  !IM
     128  !ym      INTEGER klevm1
     129  !---------------------------------------------------------------
     130  !
     131  ! Fonctions en ligne:
     132  !
     133  REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
     134  REAL zzz
     135  include "YOETHF.h"
     136  include "FCTTRE.h"
     137  fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con
     138  fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc
     139  !
     140  DATA appel1er /.TRUE./
     141  !ym
     142  zdelq=0.0
     143
     144  print*,'NUAGES4 A. JAM'
     145  IF (appel1er) THEN
     146     !
     147     PRINT*, 'fisrtilp, ninter:', ninter
     148     PRINT*, 'fisrtilp, evap_prec:', evap_prec
     149     PRINT*, 'fisrtilp, cpartiel:', cpartiel
     150     IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
     151        PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
     152        PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
     153        !         CALL abort
     154     ENDIF
     155     appel1er = .FALSE.
     156     !
     157     !AA initialiation provisoire
     158     a_tr_sca(1) = -0.5
     159     a_tr_sca(2) = -0.5
     160     a_tr_sca(3) = -0.5
     161     a_tr_sca(4) = -0.5
     162     !
     163     !AA Initialisation a 1 des coefs des fractions lessivees
     164     !
     165     !cdir collapse
     166     DO k = 1, klev
     167        DO i = 1, klon
     168           pfrac_nucl(i,k)=1.
     169           pfrac_1nucl(i,k)=1.
     170           pfrac_impa(i,k)=1.
     171        ENDDO
     172     ENDDO
     173
     174  ENDIF          !  test sur appel1er
     175  !
     176  !MAf Initialisation a 0 de zoliq
     177  !      DO i = 1, klon
     178  !         zoliq(i)=0.
     179  !      ENDDO
     180  ! Determiner les nuages froids par leur temperature
     181  !  nexpo regle la raideur de la transition eau liquide / eau glace.
     182  !
     183  ztglace = RTT - 15.0
     184  nexpo = 6
     185  !cc      nexpo = 1
     186  !
     187  ! Initialiser les sorties:
     188  !
     189  !cdir collapse
     190  DO k = 1, klev+1
     191     DO i = 1, klon
     192        prfl(i,k) = 0.0
     193        psfl(i,k) = 0.0
     194     ENDDO
     195  ENDDO
     196
     197  !cdir collapse
     198  DO k = 1, klev
     199     DO i = 1, klon
     200        d_t(i,k) = 0.0
     201        d_q(i,k) = 0.0
     202        d_ql(i,k) = 0.0
     203        rneb(i,k) = 0.0
     204        radliq(i,k) = 0.0
     205        frac_nucl(i,k) = 1.
     206        frac_impa(i,k) = 1.
     207     ENDDO
     208  ENDDO
     209  DO i = 1, klon
     210     rain(i) = 0.0
     211     snow(i) = 0.0
     212     zoliq(i)=0.
     213     !     ENDDO
     214     !
     215     ! Initialiser le flux de precipitation a zero
     216     !
     217     !     DO i = 1, klon
     218     zrfl(i) = 0.0
     219     zneb(i) = seuil_neb
     220  ENDDO
     221  !
     222  !
     223  !AA Pour plus de securite
     224
     225  zalpha_tr   = 0.
     226  zfrac_lessi = 0.
     227
     228  !AA----------------------------------------------------------
     229  !
     230  ncoreczq=0
     231  ! Boucle verticale (du haut vers le bas)
     232  !
     233  !IM : klevm1
     234  !ym      klevm1=klev-1
     235  DO k = klev, 1, -1
     236     !
     237     !AA----------------------------------------------------------
     238     !
     239     DO i = 1, klon
     240        zt(i)=t(i,k)
     241        zq(i)=q(i,k)
     242     ENDDO
     243     !
     244     ! Calculer la varition de temp. de l'air du a la chaleur sensible
     245     ! transporter par la pluie.
     246     ! Il resterait a rajouter cet effet de la chaleur sensible sur les
     247     ! flux de surface, du a la diff. de temp. entre le 1er niveau et la
     248     ! surface.
     249     !
     250     IF(k.LE.klevm1) THEN         
     251        DO i = 1, klon
     252           !IM
     253           zmair=(paprs(i,k)-paprs(i,k+1))/RG
     254           zcpair=RCPD*(1.0+RVTMP2*zq(i))
     255           zcpeau=RCPD*RVTMP2
     256           zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau &
     257                + zmair*zcpair*zt(i) ) &
     258                / (zmair*zcpair + zrfl(i)*dtime*zcpeau)
     259           !     C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
     260        ENDDO
     261     ENDIF
     262     !
     263     !
     264     ! Calculer l'evaporation de la precipitation
     265     !
     266
     267
     268     IF (evap_prec) THEN
     269        DO i = 1, klon
     270           IF (zrfl(i) .GT.0.) THEN
     271              IF (thermcep) THEN
     272                 zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
     273                 zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
     274                 zqs(i)=MIN(0.5,zqs(i))
     275                 zcor=1./(1.-RETV*zqs(i))
     276                 zqs(i)=zqs(i)*zcor
     277              ELSE
     278                 IF (zt(i) .LT. t_coup) THEN
     279                    zqs(i) = qsats(zt(i)) / pplay(i,k)
     280                 ELSE
     281                    zqs(i) = qsatl(zt(i)) / pplay(i,k)
     282                 ENDIF
     283              ENDIF
     284              zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
     285              zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
     286                   * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
     287              zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
     288                   * RG*dtime/(paprs(i,k)-paprs(i,k+1))
     289              zqev = MIN (zqev, zqevt)
     290              zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
     291                   /RG/dtime
     292
     293              ! pour la glace, on ré-évapore toute la précip dans la
     294              ! couche du dessous
     295              ! la glace venant de la couche du dessus est simplement
     296              ! dans la couche du dessous.
     297
     298              IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
     299
     300              zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
     301                   * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
     302              zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
     303                   * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
     304                   * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
     305              zrfl(i) = zrfln(i)
    282306           ENDIF
    283          ENDIF
    284          zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
    285          zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i))
    286      .         * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
    287          zqevt = MAX(0.0,MIN(zqevt,zrfl(i)))
    288      .         * RG*dtime/(paprs(i,k)-paprs(i,k+1))
    289          zqev = MIN (zqev, zqevt)
    290          zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))
    291      .                            /RG/dtime
    292 
    293 c pour la glace, on r�vapore toute la pr�ip dans la couche du dessous
    294 c la glace venant de la couche du dessus est simplement dans la couche
    295 c du dessous.
    296 
    297          IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
    298 
    299          zq(i) = zq(i) - (zrfln(i)-zrfl(i))
    300      .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
    301          zt(i) = zt(i) + (zrfln(i)-zrfl(i))
    302      .             * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
    303      .             * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
    304          zrfl(i) = zrfln(i)
    305       ENDIF
    306       ENDDO
    307       ENDIF
    308 c
    309 c Calculer Qs et L/Cp*dQs/dT:
    310 c
    311       IF (thermcep) THEN
    312          DO i = 1, klon
     307        ENDDO
     308     ENDIF
     309     !
     310     ! Calculer Qs et L/Cp*dQs/dT:
     311     !
     312     IF (thermcep) THEN
     313        DO i = 1, klon
    313314           zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
    314315           zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta
    … …  
    319320           zqs(i) = zqs(i)*zcor
    320321           zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
    321          ENDDO
    322       ELSE
    323          DO i = 1, klon
    324             IF (zt(i).LT.t_coup) THEN
    325                zqs(i) = qsats(zt(i))/pplay(i,k)
    326                zdqs(i) = dqsats(zt(i),zqs(i))
    327             ELSE
    328                zqs(i) = qsatl(zt(i))/pplay(i,k)
    329                zdqs(i) = dqsatl(zt(i),zqs(i))
    330             ENDIF
    331          ENDDO
    332       ENDIF
    333 c
    334 c Determiner la condensation partielle et calculer la quantite
    335 c de l'eau condensee:
    336 c
    337 
    338       IF (cpartiel) THEN
    339 
    340 c        print*,'Dans partiel k=',k
    341 c
    342 c   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
    343 c   nuageuse a partir des PDF de Sandrine Bony.
    344 c   rneb  : fraction nuageuse
    345 c   zqn   : eau totale dans le nuage
    346 c   zcond : eau condensee moyenne dans la maille.
    347 c           on prend en compte le r�hauffement qui diminue la partie condensee
    348 c
    349 c   Version avec les raqts
    350 
    351          if (iflag_pdf.eq.0) then
     322        ENDDO
     323     ELSE
     324        DO i = 1, klon
     325           IF (zt(i).LT.t_coup) THEN
     326              zqs(i) = qsats(zt(i))/pplay(i,k)
     327              zdqs(i) = dqsats(zt(i),zqs(i))
     328           ELSE
     329              zqs(i) = qsatl(zt(i))/pplay(i,k)
     330              zdqs(i) = dqsatl(zt(i),zqs(i))
     331           ENDIF
     332        ENDDO
     333     ENDIF
     334     !
     335     ! Determiner la condensation partielle et calculer la quantite
     336     ! de l'eau condensee:
     337     !
     338
     339     IF (cpartiel) THEN
     340
     341        !        print*,'Dans partiel k=',k
     342        !
     343        !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
     344        !   nuageuse a partir des PDF de Sandrine Bony.
     345        !   rneb  : fraction nuageuse
     346        !   zqn   : eau totale dans le nuage
     347        !   zcond : eau condensee moyenne dans la maille.
     348        !  on prend en compte le réchauffement qui diminue la partie
     349        ! condensee
     350        !
     351        !   Version avec les raqts
     352
     353        if (iflag_pdf.eq.0) then
    352354
    353355           do i=1,klon
    354             zdelq = min(ratqs(i,k),0.99) * zq(i)
    355             rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
    356             zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
     356              zdelq = min(ratqs(i,k),0.99) * zq(i)
     357              rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
     358              zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
    357359           enddo
    358360
    359          else
    360 c
    361 c   Version avec les nouvelles PDFs.
     361        else
     362           !
     363           !   Version avec les nouvelles PDFs.
    362364           do i=1,klon
    363365              if(zq(i).lt.1.e-15) then
    364                 ncoreczq=ncoreczq+1
    365                 zq(i)=1.e-15
     366                 ncoreczq=ncoreczq+1
     367                 zq(i)=1.e-15
    366368              endif
    367            enddo 
    368 
    369               if (iflag_cldcon>=5) then
    370 
    371                  call cloudth(klon,klev,k,ztv,
    372      .           zq,zqta,fraca,
    373      .           qcloud,ctot,zpspsk,paprs,ztla,zthl,
    374      .           ratqs,zqs,t)
    375 
    376                  do i=1,klon
     369           enddo
     370
     371           if (iflag_cldcon>=5) then
     372
     373              call cloudth(klon,klev,k,ztv, &
     374                   zq,zqta,fraca, &
     375                   qcloud,ctot,zpspsk,paprs,ztla,zthl, &
     376                   ratqs,zqs,t)
     377
     378              do i=1,klon
    377379                 rneb(i,k)=ctot(i,k)
    378380                 zqn(i)=qcloud(i)
    379                  enddo
    380 
     381              enddo
     382
     383           endif
     384
     385           if (iflag_cldcon <= 4) then
     386              lognormale = .true.
     387           elseif (iflag_cldcon == 6) then
     388              ! lognormale en l'absence des thermiques
     389              lognormale = fraca(:,k) < 1e-10
     390           else
     391              ! Dans le cas iflag_cldcon=5, on prend systématiquement la
     392              ! bi-gaussienne
     393              lognormale = .false.
     394           end if
     395
     396           do i=1,klon
     397              if (lognormale(i)) then
     398                 zpdf_sig(i)=ratqs(i,k)*zq(i)
     399                 zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
     400                 zpdf_delta(i)=log(zq(i)/zqs(i))
     401                 zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
     402                 zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
     403                 zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
     404                 zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
     405                 zpdf_e1(i)=1.-erf(zpdf_e1(i))
     406                 zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
     407                 zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
     408                 zpdf_e2(i)=1.-erf(zpdf_e2(i))
    381409              endif
    382 
    383 ! Pour iflag_cldcon<=4, on prend toujours la lognormale
    384 ! Dans le cas iflag_cldcon=5, on prend systématiquement la bi-gaussienne
    385 ! Dans le cas iflagÃ_cldcon=6, on prend la lognormale en absence des thermiques
    386 
    387             lognormale(:)=
    388      .      iflag_cldcon<=4.or.(iflag_cldcon==6.and.fraca(:,k)<1.e-10)
    389             do i=1,klon
    390             if (lognormale(i)) then
    391             zpdf_sig(i)=ratqs(i,k)*zq(i)
    392             zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
    393             zpdf_delta(i)=log(zq(i)/zqs(i))
    394             zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
    395             zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
    396             zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
    397             zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
    398             zpdf_e1(i)=1.-erf(zpdf_e1(i))
    399             zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
    400             zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
    401             zpdf_e2(i)=1.-erf(zpdf_e2(i))
    402             endif
    403             enddo
    404 
    405             do i=1,klon
    406             if (lognormale(i)) then
    407             if (zpdf_e1(i).lt.1.e-10) then
    408                rneb(i,k)=0.
    409                zqn(i)=zqs(i)
    410             else
    411                rneb(i,k)=0.5*zpdf_e1(i)
    412                zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
    413             endif
    414             endif
    415            
    416            enddo
     410           enddo
     411
     412           do i=1,klon
     413              if (lognormale(i)) then
     414                 if (zpdf_e1(i).lt.1.e-10) then
     415                    rneb(i,k)=0.
     416                    zqn(i)=zqs(i)
     417                 else
     418                    rneb(i,k)=0.5*zpdf_e1(i)
     419                    zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
     420                 endif
     421              endif
     422
     423           enddo
    417424
    418425
    … …  
    435442           ENDIF
    436443        ENDDO
    437 !         do i=1,klon
    438 !            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
    439 !            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
    440 !            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
    441 !c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
    442 !c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
    443 !c  la convection.
    444 !c  ATTENTION !!! Il va falloir verifier tout ca.
    445 !            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
    446 !c           print*,'ZDQS ',zdqs(i)
    447 !c--Olivier
    448 !            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
    449 !            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
    450 !            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
    451 !c--fin
    452 !           ENDDO
    453       ELSE
    454          DO i = 1, klon
    455             IF (zq(i).GT.zqs(i)) THEN
    456                rneb(i,k) = 1.0
    457             ELSE
    458                rneb(i,k) = 0.0
    459             ENDIF
    460             zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
    461          ENDDO
    462       ENDIF
    463 c
    464       DO i = 1, klon
    465          zq(i) = zq(i) - zcond(i)
    466 c         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
    467          zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
    468       ENDDO
    469 c
    470 c Partager l'eau condensee en precipitation et eau liquide nuageuse
    471 c
    472       DO i = 1, klon
    473       IF (rneb(i,k).GT.0.0) THEN
    474          zoliq(i) = zcond(i)
    475          zrho(i) = pplay(i,k) / zt(i) / RD
    476          zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
    477          zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
    478          zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
    479          zfice(i) = zfice(i)**nexpo
    480          zneb(i) = MAX(rneb(i,k), seuil_neb)
    481          radliq(i,k) = zoliq(i)/REAL(ninter+1)
    482       ENDIF
    483       ENDDO
    484 c
    485       DO n = 1, ninter
    486       DO i = 1, klon
    487       IF (rneb(i,k).GT.0.0) THEN
    488          zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
    489 
    490          IF (zneb(i).EQ.seuil_neb) THEN
    491              ztot = 0.0
    492          ELSE
    493 c  quantite d'eau a eliminer: zchau
    494 c  meme chose pour la glace: zfroi
    495              if (ptconv(i,k)) then
    496                 zcl   =cld_lc_con
    497                 zct   =1./cld_tau_con
    498                 zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i)
    499      .              *fallvc(zrhol(i)) * zfice(i)
    500              else
    501                 zcl   =cld_lc_lsc
    502                 zct   =1./cld_tau_lsc
    503                 zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i)
    504      .              *fallvs(zrhol(i)) * zfice(i)
    505              endif
    506              zchau    = zct   *dtime/REAL(ninter) * zoliq(i)
    507      .         *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
    508              ztot    = zchau    + zfroi
    509              ztot    = MAX(ztot   ,0.0)
    510          ENDIF
    511          ztot    = MIN(ztot,zoliq(i))
    512          zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
    513          radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
    514       ENDIF
    515       ENDDO
    516       ENDDO
    517 c
    518       DO i = 1, klon
    519       IF (rneb(i,k).GT.0.0) THEN
    520          d_ql(i,k) = zoliq(i)
    521          zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0)
    522      .                    * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
    523       ENDIF
    524       IF (zt(i).LT.RTT) THEN
    525         psfl(i,k)=zrfl(i)
    526       ELSE
    527         prfl(i,k)=zrfl(i)
    528       ENDIF
    529       ENDDO
    530 c
    531 c Calculer les tendances de q et de t:
    532 c
    533       DO i = 1, klon
    534          d_q(i,k) = zq(i) - q(i,k)
    535          d_t(i,k) = zt(i) - t(i,k)
    536       ENDDO
    537 c
    538 cAA--------------- Calcul du lessivage stratiforme  -------------
    539 
    540       DO i = 1,klon
    541 c
    542          zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)
    543      .                * (paprs(i,k)-paprs(i,k+1))/RG
    544          IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
    545 cAA lessivage nucleation LMD5 dans la couche elle-meme
    546             if (t(i,k) .GE. ztglace) THEN
    547                zalpha_tr = a_tr_sca(3)
    548             else
    549                zalpha_tr = a_tr_sca(4)
    550             endif
    551             zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
    552             pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
    553             frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
    554 c
    555 c nucleation avec un facteur -1 au lieu de -0.5
    556             zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
    557             pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
    558          ENDIF
    559 c
    560       ENDDO      ! boucle sur i
    561 c
    562 cAA Lessivage par impaction dans les couches en-dessous
    563       DO kk = k-1, 1, -1
    564         DO i = 1, klon
    565           IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
    566             if (t(i,kk) .GE. ztglace) THEN
    567               zalpha_tr = a_tr_sca(1)
    568             else
    569               zalpha_tr = a_tr_sca(2)
    570             endif
    571             zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
    572             pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
    573             frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
    574           ENDIF
    575         ENDDO
    576       ENDDO
    577 c
    578 cAA----------------------------------------------------------
    579 c                     FIN DE BOUCLE SUR K   
    580  9999 CONTINUE
    581 c
    582 cAA-----------------------------------------------------------
    583 c
    584 c Pluie ou neige au sol selon la temperature de la 1ere couche
    585 c
    586       DO i = 1, klon
    587       IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
    588          snow(i) = zrfl(i)
    589          zlh_solid(i) = RLSTT-RLVTT
    590       ELSE
    591          rain(i) = zrfl(i)
    592          zlh_solid(i) = 0.
    593       ENDIF
    594       ENDDO
    595 C
    596 C For energy conservation : when snow is present, the solification
    597 c latent heat is considered.
    598       DO k = 1, klev
    599         DO i = 1, klon
    600           zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
    601           zmair=(paprs(i,k)-paprs(i,k+1))/RG
    602           zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
    603           d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
    604         END DO
    605       END DO
    606 c
    607 
    608       if (ncoreczq>0) then
    609          print*,'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
    610       endif
    611       RETURN
    612       END
     444        !         do i=1,klon
     445        !            IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
     446        !            IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
     447        !            rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
     448        !c           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
     449        !c  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
     450        !c  la convection.
     451        !c  ATTENTION !!! Il va falloir verifier tout ca.
     452        !            zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
     453        !c           print*,'ZDQS ',zdqs(i)
     454        !c--Olivier
     455        !            rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
     456        !            IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
     457        !            IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
     458        !c--fin
     459        !           ENDDO
     460     ELSE
     461        DO i = 1, klon
     462           IF (zq(i).GT.zqs(i)) THEN
     463              rneb(i,k) = 1.0
     464           ELSE
     465              rneb(i,k) = 0.0
     466           ENDIF
     467           zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
     468        ENDDO
     469     ENDIF
     470     !
     471     DO i = 1, klon
     472        zq(i) = zq(i) - zcond(i)
     473        !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
     474        zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
     475     ENDDO
     476     !
     477     ! Partager l'eau condensee en precipitation et eau liquide nuageuse
     478     !
     479     DO i = 1, klon
     480        IF (rneb(i,k).GT.0.0) THEN
     481           zoliq(i) = zcond(i)
     482           zrho(i) = pplay(i,k) / zt(i) / RD
     483           zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
     484           zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
     485           zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
     486           zfice(i) = zfice(i)**nexpo
     487           zneb(i) = MAX(rneb(i,k), seuil_neb)
     488           radliq(i,k) = zoliq(i)/REAL(ninter+1)
     489        ENDIF
     490     ENDDO
     491     !
     492     DO n = 1, ninter
     493        DO i = 1, klon
     494           IF (rneb(i,k).GT.0.0) THEN
     495              zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
     496
     497              IF (zneb(i).EQ.seuil_neb) THEN
     498                 ztot = 0.0
     499              ELSE
     500                 !  quantite d'eau a eliminer: zchau
     501                 !  meme chose pour la glace: zfroi
     502                 if (ptconv(i,k)) then
     503                    zcl   =cld_lc_con
     504                    zct   =1./cld_tau_con
     505                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
     506                         *fallvc(zrhol(i)) * zfice(i)
     507                 else
     508                    zcl   =cld_lc_lsc
     509                    zct   =1./cld_tau_lsc
     510                    zfroi    = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
     511                         *fallvs(zrhol(i)) * zfice(i)
     512                 endif
     513                 zchau    = zct   *dtime/REAL(ninter) * zoliq(i) &
     514                      *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl   )**2)) *(1.-zfice(i))
     515                 ztot    = zchau    + zfroi
     516                 ztot    = MAX(ztot   ,0.0)
     517              ENDIF
     518              ztot    = MIN(ztot,zoliq(i))
     519              zoliq(i) = MAX(zoliq(i)-ztot   , 0.0)
     520              radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
     521           ENDIF
     522        ENDDO
     523     ENDDO
     524     !
     525     DO i = 1, klon
     526        IF (rneb(i,k).GT.0.0) THEN
     527           d_ql(i,k) = zoliq(i)
     528           zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
     529                * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
     530        ENDIF
     531        IF (zt(i).LT.RTT) THEN
     532           psfl(i,k)=zrfl(i)
     533        ELSE
     534           prfl(i,k)=zrfl(i)
     535        ENDIF
     536     ENDDO
     537     !
     538     ! Calculer les tendances de q et de t:
     539     !
     540     DO i = 1, klon
     541        d_q(i,k) = zq(i) - q(i,k)
     542        d_t(i,k) = zt(i) - t(i,k)
     543     ENDDO
     544     !
     545     !AA--------------- Calcul du lessivage stratiforme  -------------
     546
     547     DO i = 1,klon
     548        !
     549        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
     550             * (paprs(i,k)-paprs(i,k+1))/RG
     551        IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
     552           !AA lessivage nucleation LMD5 dans la couche elle-meme
     553           if (t(i,k) .GE. ztglace) THEN
     554              zalpha_tr = a_tr_sca(3)
     555           else
     556              zalpha_tr = a_tr_sca(4)
     557           endif
     558           zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
     559           pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
     560           frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
     561           !
     562           ! nucleation avec un facteur -1 au lieu de -0.5
     563           zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
     564           pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)
     565        ENDIF
     566        !
     567     ENDDO      ! boucle sur i
     568     !
     569     !AA Lessivage par impaction dans les couches en-dessous
     570     DO kk = k-1, 1, -1
     571        DO i = 1, klon
     572           IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
     573              if (t(i,kk) .GE. ztglace) THEN
     574                 zalpha_tr = a_tr_sca(1)
     575              else
     576                 zalpha_tr = a_tr_sca(2)
     577              endif
     578              zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
     579              pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)
     580              frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
     581           ENDIF
     582        ENDDO
     583     ENDDO
     584     !
     585     !AA----------------------------------------------------------
     586     !                     FIN DE BOUCLE SUR K   
     587  end DO
     588  !
     589  !AA-----------------------------------------------------------
     590  !
     591  ! Pluie ou neige au sol selon la temperature de la 1ere couche
     592  !
     593  DO i = 1, klon
     594     IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
     595        snow(i) = zrfl(i)
     596        zlh_solid(i) = RLSTT-RLVTT
     597     ELSE
     598        rain(i) = zrfl(i)
     599        zlh_solid(i) = 0.
     600     ENDIF
     601  ENDDO
     602  !
     603  ! For energy conservation : when snow is present, the solification
     604  ! latent heat is considered.
     605  DO k = 1, klev
     606     DO i = 1, klon
     607        zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k)))
     608        zmair=(paprs(i,k)-paprs(i,k+1))/RG
     609        zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
     610        d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair)
     611     END DO
     612  END DO
     613  !
     614
     615  if (ncoreczq>0) then
     616     print*,'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
     617  endif
     618
     619END SUBROUTINE fisrtilp
Note: See TracChangeset for help on using the changeset viewer.