Ignore:
Timestamp:
Oct 14, 2016, 2:57:28 PM (8 years ago)
Author:
Laurent Fairhead
Message:

Merged trunk changes r2640:2664 into testing branch

Location:
LMDZ5/branches/testing
Files:
1 deleted
26 edited
2 copied

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/testing

  • LMDZ5/branches/testing/libf/dynphy_lonlat/phydev/callphysiq_mod.F90

    r2435 r2669  
    1212                       jD_cur,jH_cur_split,zdt_split,                     &
    1313                       zplev_omp,zplay_omp,                               &
    14                        zphi_omp,zphis_omp,                                &
     14                       zpk_omp,zphi_omp,zphis_omp,                        &
    1515                       presnivs_omp,                                      &
    1616                       zufi_omp,zvfi_omp,zrfi_omp,ztfi_omp,zqfi_omp,      &
    … …  
    3434  REAL,INTENT(IN) :: zplev_omp(klon,llm+1) ! interlayer pressure (Pa)
    3535  REAL,INTENT(IN) :: zplay_omp(klon,llm) ! mid-layer pressure (Pa)
     36  REAL,INTENT(IN) :: zpk_omp(klon,llm) ! Exner function
    3637  REAL,INTENT(IN) :: zphi_omp(klon,llm) ! geopotential at midlayer
    3738  REAL,INTENT(IN) :: zphis_omp(klon) ! surface geopotential

  • LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90

    r2641 r2669  
    111111  INTEGER :: flag_aerosol
    112112  INTEGER :: flag_aerosol_strat
     113  LOGICAL :: flag_bc_internal_mixture
    113114  LOGICAL :: new_aod
    114115  REAL    :: bl95_b0, bl95_b1
    … …  
    131132                   ok_ade, ok_aie, ok_cdnc, aerosol_couple,             &
    132133                   flag_aerosol, flag_aerosol_strat, new_aod,           &
    133                    bl95_b0, bl95_b1,                                    &
     134                   flag_bc_internal_mixture, bl95_b0, bl95_b1,          &
    134135                   read_climoz,                                         &
    135136                   alp_offset)

  • LMDZ5/branches/testing/libf/dynphy_lonlat/phylmd/iniphysiq_mod.F90

    r2641 r2669  
    99                     nbp, communicator, &
    1010                     punjours, pdayref,ptimestep, &
    11                      rlatu,rlatv,rlonu,rlonv,aire,cu,cv,      &
     11                     rlatudyn,rlatvdyn,rlonudyn,rlonvdyn,airedyn,cudyn,cvdyn, &
    1212                     prad,pg,pr,pcpp,iflag_phys)
    1313  USE dimphy, ONLY: init_dimphy
    … …  
    4444  USE mod_phys_lmdz_omp_data, ONLY: klon_omp
    4545#endif
     46  USE ioipsl_getin_p_mod, ONLY: getin_p
     47  USE slab_heat_transp_mod, ONLY: ini_slab_transp_geom
    4648  IMPLICIT NONE
    4749
    … …  
    5254
    5355  include "dimensions.h"
     56  include "paramet.h"
    5457  include "iniprint.h"
    5558  include "tracstoke.h"
     59  include "comgeom.h"
    5660
    5761  REAL, INTENT (IN) :: prad ! radius of the planet (m)
    … …  
    6569  INTEGER, INTENT(IN) :: nbp ! number of physics columns for this MPI process
    6670  INTEGER, INTENT(IN) :: communicator ! MPI communicator
    67   REAL, INTENT (IN) :: rlatu(jj+1) ! latitudes of the physics grid
    68   REAL, INTENT (IN) :: rlatv(jj) ! latitude boundaries of the physics grid
    69   REAL, INTENT (IN) :: rlonv(ii+1) ! longitudes of the physics grid
    70   REAL, INTENT (IN) :: rlonu(ii+1) ! longitude boundaries of the physics grid
    71   REAL, INTENT (IN) :: aire(ii+1,jj+1) ! area of the dynamics grid (m2)
    72   REAL, INTENT (IN) :: cu((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u)
    73   REAL, INTENT (IN) :: cv((ii+1)*jj) ! cv coeff. (v_covariant = cv * v)
     71  REAL, INTENT (IN) :: rlatudyn(jj+1) ! latitudes of the physics grid
     72  REAL, INTENT (IN) :: rlatvdyn(jj) ! latitude boundaries of the physics grid
     73  REAL, INTENT (IN) :: rlonvdyn(ii+1) ! longitudes of the physics grid
     74  REAL, INTENT (IN) :: rlonudyn(ii+1) ! longitude boundaries of the physics grid
     75  REAL, INTENT (IN) :: airedyn(ii+1,jj+1) ! area of the dynamics grid (m2)
     76  REAL, INTENT (IN) :: cudyn((ii+1)*(jj+1)) ! cu coeff. (u_covariant = cu * u)
     77  REAL, INTENT (IN) :: cvdyn((ii+1)*jj) ! cv coeff. (v_covariant = cv * v)
    7478  INTEGER, INTENT (IN) :: pdayref ! reference day of for the simulation
    7579  REAL, INTENT (IN) :: ptimestep !physics time step (s)
    … …  
    8084  CHARACTER (LEN=20) :: modname = 'iniphysiq'
    8185  CHARACTER (LEN=80) :: abort_message
    82 
     86 
     87  LOGICAL :: slab_hdiff
     88  INTEGER :: slab_ekman
     89  CHARACTER (LEN = 6) :: type_ocean
    8390
    8491#ifndef CPP_PARA
    … …  
    9299  CALL inigeomphy(ii,jj,nlayer, &
    93100               nbp, communicator, &
    94                rlatu,rlatv, &
    95                rlonu,rlonv, &
    96                aire,cu,cv)
     101               rlatudyn,rlatvdyn, &
     102               rlonudyn,rlonvdyn, &
     103               airedyn,cudyn,cvdyn)
    97104
    98105  ! --> now initialize things specific to the phylmd physics package
    … …  
    117124  ! Copy over "offline" settings
    118125  CALL init_phystokenc(offline,istphy)
     126
     127  ! Initialization for slab heat transport
     128  type_ocean="force"
     129  CALL getin_p('type_ocean',type_ocean)
     130  slab_hdiff=.FALSE.
     131  CALL getin_p('slab_hdiff',slab_hdiff)
     132  slab_ekman=0
     133  CALL getin_p('slab_ekman',slab_ekman)
     134  IF ((type_ocean=='slab').AND.(slab_hdiff.OR.(slab_ekman.GT.0))) THEN
     135     CALL ini_slab_transp_geom(ip1jm,ip1jmp1,unsairez,fext,unsaire,&
     136                                  cu,cuvsurcv,cv,cvusurcu, &
     137                                  aire,apoln,apols, &
     138                                  aireu,airev)
     139  END IF
    119140
    120141  ! Initialize tracer names, numbers, etc. for physics
    … …  
    159180#ifdef INCA
    160181     CALL init_inca_dim(klon_omp,nbp_lev,nbp_lon,nbp_lat - 1, &
    161           rlonu,rlatu,rlonv,rlatv)
     182          rlonudyn,rlatudyn,rlonvdyn,rlatvdyn)
    162183#endif
    163184  END IF

  • LMDZ5/branches/testing/libf/phydev/physiq_mod.F90

    r2435 r2669  
    2525#ifdef CPP_XIOS
    2626      USE xios, ONLY: xios_update_calendar
    27       USE wxios, only: wxios_add_vaxis, wxios_set_timestep, wxios_closedef, &
    28                       histwrite_phy
     27      USE wxios, only: wxios_add_vaxis, wxios_set_cal, wxios_closedef
     28      USE iophy, ONLY: histwrite_phy
    2929#endif
    3030
    … …  
    129129!XIOS
    130130    ! Declare available vertical axes to be used in output files:   
    131     !CALL wxios_add_vaxis("presnivs", "dummy-not-used", klev, presnivs)
    132131    CALL wxios_add_vaxis("presnivs", klev, presnivs)
    133132
    134     ! Declare time step length (in s):
    135     CALL wxios_set_timestep(dtime)
    136 
     133    ! Declare calendar and time step
     134    CALL wxios_set_cal(dtime,"earth_360d",1,1,1,0.0,1,1,1,0.0)
     135   
    137136    !Finalize the context:
    138137    CALL wxios_closedef()
    139138#endif
    140139!$OMP END MASTER
     140!$OMP BARRIER
    141141endif ! of if (debut)
    142142

  • LMDZ5/branches/testing/libf/phylmd/Dust/phytracr_spl_mod.F90

    r2641 r2669  
    11701170       !print *,nbtr
    11711171       do it=1,nbtr
    1172         print *, it, tname(it+2)
    1173         if (tname(it+2) == 'PREC' ) then
     1172        print *, it, tname(it+nqo)
     1173        if (tname(it+nqo) == 'PREC' ) then
    11741174            id_prec=it
    11751175        endif
    1176         if (tname(it+2) == 'FINE' ) then
     1176        if (tname(it+nqo) == 'FINE' ) then
    11771177            id_fine=it
    11781178        endif
    1179         if (tname(it+2) == 'COSS' ) then
     1179        if (tname(it+nqo) == 'COSS' ) then
    11801180            id_coss=it
    11811181        endif
    1182         if (tname(it+2) == 'CODU' ) then
     1182        if (tname(it+nqo) == 'CODU' ) then
    11831183            id_codu=it
    11841184        endif
    1185         if (tname(it+2) == 'SCDU' ) then
     1185        if (tname(it+nqo) == 'SCDU' ) then
    11861186            id_scdu=it
    11871187        endif
    … …  
    26212621         call iophys_ecrit('fav'//str2,1,'SOURCE','',source_tr(:,it))
    26222622      enddo
    2623 #endif
    2624 
    26252623      do it=1,nbtr
    26262624         write(str2,'(i2.2)') it
    26272625         call iophys_ecrit('TRB'//str2,klev,'SOURCE','',tr_seri(:,:,it))
    26282626      enddo
     2627#endif
    26292628
    26302629
    … …  
    26502649
    26512650#ifdef IOPHYS_DUST
    2652 
    26532651      do it=1,nbtr
    26542652         write(str2,'(i2.2)') it
    … …  
    26762674!
    26772675#ifdef IOPHYS_DUST
    2678 !
    26792676      print *,'INPUT TO PRECUREMISSION'
    26802677         call iophys_ecrit('ftsol',4,'ftsol','',ftsol)

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

    r2408 r2669  
    8787       CALL limit_read_frac(itime, dtime, jour, pctsrf, is_modified)
    8888    CASE ('slab')
     89       IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN
     90       ! Read fraction from limit.nc
     91           CALL limit_read_frac(itime, dtime, jour, pctsrf, is_modified)
     92       ELSE
    8993       ! Get fraction from slab module
    90        CALL ocean_slab_frac(itime, dtime, jour, pctsrf, is_modified)
     94           CALL ocean_slab_frac(itime, dtime, jour, pctsrf, is_modified)
     95       ENDIF
    9196    CASE ('couple')
    9297       ! Get fraction from the coupler

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

    r2594 r2669  
    99
    1010!===========================================================================
    11 ! Auteur : Arnaud Octavio Jam (LMD/CNRS)
     11! Author : Arnaud Octavio Jam (LMD/CNRS)
    1212! Date : 25 Mai 2010
    1313! Objet : calcule les valeurs de qc et rneb dans les thermiques
    … …  
    4949      REAL rneb(ngrid,klev)
    5050      REAL t(ngrid,klev)
    51       REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
     51      REAL qsatmmussig1,qsatmmussig2,sqrt2pi,sqrt2,sqrtpi,pi
    5252      REAL rdd,cppd,Lv
    5353      REAL alth,alenv,ath,aenv
    54       REAL sth,senv,sigma1s,sigma2s,xth,xenv
     54      REAL sth,senv,sigma1s,sigma2s,xth,xenv, exp_xenv1, exp_xenv2,exp_xth1,exp_xth2
    5555      REAL Tbef,zdelta,qsatbef,zcor
    5656      REAL qlbef 
    57       REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
    58      
     57      REAL ratqs(ngrid,klev) ! Determine the width of the vapour distribution
    5958      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
    6059      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
    … …  
    6362
    6463
    65 
    66 
    67 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    68 ! Gestion de deux versions de cloudth
    69 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    7064
    7165      IF (iflag_cloudth_vert.GE.1) THEN
    … …  
    8276! Initialisation des variables r?elles
    8377!-------------------------------------------------------------------------------
    84       sigma1(:,:)=0.
    85       sigma2(:,:)=0.
    86       qlth(:,:)=0.
    87       qlenv(:,:)=0. 
    88       qltot(:,:)=0.
    89       rneb(:,:)=0.
    90       qcloud(:)=0.
    91       cth(:,:)=0.
    92       cenv(:,:)=0.
    93       ctot(:,:)=0.
    94       qsatmmussig1=0.
    95       qsatmmussig2=0.
    96       rdd=287.04
    97       cppd=1005.7
    98       pi=3.14159
    99       Lv=2.5e6
    100       sqrt2pi=sqrt(2.*pi)
    101 
    102 
    103 
    104 !-------------------------------------------------------------------------------
    105 ! Calcul de la fraction du thermique et des ?cart-types des distributions
    106 !-------------------------------------------------------------------------------                 
    107       do ind1=1,ngrid
    108 
    109       if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then
    110 
    111       zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
    112 
    113 
    114 !      zqenv(ind1)=po(ind1)
    115       Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
    116       zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    117       qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
    118       qsatbef=MIN(0.5,qsatbef)
    119       zcor=1./(1.-retv*qsatbef)
    120       qsatbef=qsatbef*zcor
    121       zqsatenv(ind1,ind2)=qsatbef
    122 
    123 
    124 
    125 
    126       alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
    127       aenv=1./(1.+(alenv*Lv/cppd))
    128       senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
    129 
    130 
    131 
    132 
    133       Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
    134       zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    135       qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
    136       qsatbef=MIN(0.5,qsatbef)
    137       zcor=1./(1.-retv*qsatbef)
    138       qsatbef=qsatbef*zcor
    139       zqsatth(ind1,ind2)=qsatbef
    140            
    141       alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
    142       ath=1./(1.+(alth*Lv/cppd))
    143       sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
    144      
    145      
    146 
    147 !------------------------------------------------------------------------------
    148 ! Calcul des ?cart-types pour s
    149 !------------------------------------------------------------------------------
    150 
    151 !      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
    152 !      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.002*zqta(ind1,ind2)
    153 !       if (paprs(ind1,ind2).gt.90000) then
    154 !       ratqs(ind1,ind2)=0.002
    155 !       else
    156 !       ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
    157 !       endif
    158        sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
    159        sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
    160 !       sigma1s=ratqs(ind1,ind2)*po(ind1)
    161 !      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
    162  
    163 !------------------------------------------------------------------------------
    164 ! Calcul de l'eau condens?e et de la couverture nuageuse
    165 !------------------------------------------------------------------------------
    166       sqrt2pi=sqrt(2.*pi)
    167       xth=sth/(sqrt(2.)*sigma2s)
    168       xenv=senv/(sqrt(2.)*sigma1s)
    169       cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
    170       cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
    171       ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
    172 
    173       qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
    174       qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
    175       qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
    176 
    177 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    178       if (ctot(ind1,ind2).lt.1.e-10) then
    179       ctot(ind1,ind2)=0.
    180       qcloud(ind1)=zqsatenv(ind1,ind2)
    181 
    182       else   
    183                
    184       ctot(ind1,ind2)=ctot(ind1,ind2)
    185       qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
    186 
    187       endif                           
    188      
    189          
    190 !     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
    191 
    192 
    193       else  ! gaussienne environnement seule
    194      
    195       zqenv(ind1)=po(ind1)
    196       Tbef=t(ind1,ind2)
    197       zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    198       qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
    199       qsatbef=MIN(0.5,qsatbef)
    200       zcor=1./(1.-retv*qsatbef)
    201       qsatbef=qsatbef*zcor
    202       zqsatenv(ind1,ind2)=qsatbef
    203      
    204 
    205 !      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
    206       zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
    207       alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
    208       aenv=1./(1.+(alenv*Lv/cppd))
    209       senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
    210      
    211 
    212       sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
    213 
    214       sqrt2pi=sqrt(2.*pi)
    215       xenv=senv/(sqrt(2.)*sigma1s)
    216       ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
    217       qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
    218      
    219       if (ctot(ind1,ind2).lt.1.e-3) then
    220       ctot(ind1,ind2)=0.
    221       qcloud(ind1)=zqsatenv(ind1,ind2)
    222 
    223       else   
    224                
    225       ctot(ind1,ind2)=ctot(ind1,ind2)
    226       qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
    227 
    228       endif   
    229  
    230  
    231  
    232  
    233  
    234  
    235       endif   
    236       enddo
    237      
    238       return
    239       end
    240 
    241 
    242 
    243 !===========================================================================
    244      SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
    245      &           ztv,po,zqta,fraca, &
    246      &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
    247      &           ratqs,zqs,t)
    248 
    249 !===========================================================================
    250 ! Auteur : Arnaud Octavio Jam (LMD/CNRS)
    251 ! Date : 25 Mai 2010
    252 ! Objet : calcule les valeurs de qc et rneb dans les thermiques
    253 !===========================================================================
    254 
    255 
    256       USE ioipsl_getin_p_mod, ONLY : getin_p
    257 
    258       IMPLICIT NONE
    259 
    260 #include "YOMCST.h"
    261 #include "YOETHF.h"
    262 #include "FCTTRE.h"
    263 #include "thermcell.h"
    264 #include "nuage.h"
    265      
    266       INTEGER itap,ind1,ind2
    267       INTEGER ngrid,klev,klon,l,ig
    268      
    269       REAL ztv(ngrid,klev)
    270       REAL po(ngrid)
    271       REAL zqenv(ngrid)   
    272       REAL zqta(ngrid,klev)
    273          
    274       REAL fraca(ngrid,klev+1)
    275       REAL zpspsk(ngrid,klev)
    276       REAL paprs(ngrid,klev+1)
    277       REAL ztla(ngrid,klev)
    278       REAL zthl(ngrid,klev)
    279 
    280       REAL zqsatth(ngrid,klev)
    281       REAL zqsatenv(ngrid,klev)
    282      
    283      
    284       REAL sigma1(ngrid,klev)                                                         
    285       REAL sigma2(ngrid,klev)
    286       REAL qlth(ngrid,klev)
    287       REAL qlenv(ngrid,klev)
    288       REAL qltot(ngrid,klev)
    289       REAL cth(ngrid,klev) 
    290       REAL cenv(ngrid,klev)   
    291       REAL ctot(ngrid,klev)
    292       REAL rneb(ngrid,klev)
    293       REAL t(ngrid,klev)                                                                 
    294       REAL qsatmmussig1,qsatmmussig2,sqrt2pi,pi
    295       REAL rdd,cppd,Lv,sqrt2,sqrtpi
    296       REAL alth,alenv,ath,aenv
    297       REAL sth,senv,sigma1s,sigma2s,xth,xenv
    298       REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
    299       REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
    300       REAL Tbef,zdelta,qsatbef,zcor
    301       REAL qlbef 
    302       REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
    303       ! Change the width of the PDF used for vertical subgrid scale heterogeneity
    304       ! (J Jouhaud, JL Dufresne, JB Madeleine)
    305       REAL,SAVE :: vert_alpha
    306       LOGICAL, SAVE :: firstcall = .TRUE.
    307      
    308       REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
    309       REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
    310       REAL zqs(ngrid), qcloud(ngrid)
    311       REAL erf
    312 
    313 !------------------------------------------------------------------------------
    314 ! Initialisation des variables r?elles
    315 !------------------------------------------------------------------------------
    31678      sigma1(:,:)=0.
    31779      sigma2(:,:)=0.
    … …  
    33496      sqrtpi=sqrt(pi)
    33597
    336       IF (firstcall) THEN
    337         vert_alpha=0.5
    338         CALL getin_p('cloudth_vert_alpha',vert_alpha)
    339         WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
    340         firstcall=.FALSE.
    341       ENDIF
    342 
    343 !-------------------------------------------------------------------------------
    344 ! Calcul de la fraction du thermique et des ?cart-types des distributions
     98
     99!-------------------------------------------------------------------------------
     100! Cloud fraction in the thermals and standard deviation of the PDFs
    345101!-------------------------------------------------------------------------------                 
    346102      do ind1=1,ngrid
    … …  
    350106      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2))
    351107
    352 
    353 !      zqenv(ind1)=po(ind1)
    354108      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
    355109      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    356       qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
     110      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
    357111      qsatbef=MIN(0.5,qsatbef)
    358112      zcor=1./(1.-retv*qsatbef)
    … …  
    361115
    362116
    363 
    364 
    365       alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
    366       aenv=1./(1.+(alenv*Lv/cppd))
    367       senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
    368 
    369 
     117      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
     118      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
     119      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
     120
     121!po = qt de l'environnement ET des thermique
     122!zqenv = qt environnement
     123!zqsatenv = qsat environnement
     124!zthl = Tl environnement
    370125
    371126
    … …  
    378133      zqsatth(ind1,ind2)=qsatbef
    379134           
    380       alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)   
    381       ath=1./(1.+(alth*Lv/cppd))
    382       sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))
    383      
    384      
    385 
    386 !------------------------------------------------------------------------------
    387 ! Calcul des ?cart-types pour s
    388 !------------------------------------------------------------------------------
    389 
    390       sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
    391       sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
     135      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
     136      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
     137      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
     138     
     139!zqta = qt thermals
     140!zqsatth = qsat thermals
     141!ztla = Tl thermals
     142
     143!------------------------------------------------------------------------------
     144! s standard deviations
     145!------------------------------------------------------------------------------
     146
     147!     tests
     148!     sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+0.002*po(ind1)
     149!     sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+ratqs(ind1,ind2)*po(ind1)
     150!     sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
     151!     final option
     152      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
     153      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.01)**0.4+0.002*zqta(ind1,ind2)
     154 
     155!------------------------------------------------------------------------------
     156! Condensed water and cloud cover
     157!------------------------------------------------------------------------------
     158      xth=sth/(sqrt2*sigma2s)
     159      xenv=senv/(sqrt2*sigma1s)
     160      cth(ind1,ind2)=0.5*(1.+1.*erf(xth))       !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
     161      cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))     !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
     162      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
     163
     164      qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt2*cth(ind1,ind2))
     165      qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
     166      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
     167
     168      if (ctot(ind1,ind2).lt.1.e-10) then
     169      ctot(ind1,ind2)=0.
     170      qcloud(ind1)=zqsatenv(ind1,ind2)
     171      else
     172      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
     173      endif
     174
     175      else  ! Environnement only, follow the if l.110
     176     
     177      zqenv(ind1)=po(ind1)
     178      Tbef=t(ind1,ind2)
     179      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
     180      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
     181      qsatbef=MIN(0.5,qsatbef)
     182      zcor=1./(1.-retv*qsatbef)
     183      qsatbef=qsatbef*zcor
     184      zqsatenv(ind1,ind2)=qsatbef
     185
     186!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
     187      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
     188      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
     189      aenv=1./(1.+(alenv*Lv/cppd))
     190      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))
     191     
     192      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
     193
     194      xenv=senv/(sqrt2*sigma1s)
     195      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
     196      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
     197
     198      if (ctot(ind1,ind2).lt.1.e-3) then
     199      ctot(ind1,ind2)=0.
     200      qcloud(ind1)=zqsatenv(ind1,ind2)
     201      else   
     202      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
     203      endif
     204
     205
     206      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
     207      enddo       ! from the loop on ngrid l.108
     208      return
     209      end
     210
     211
     212
     213!===========================================================================
     214     SUBROUTINE cloudth_vert(ngrid,klev,ind2,  &
     215     &           ztv,po,zqta,fraca, &
     216     &           qcloud,ctot,zpspsk,paprs,ztla,zthl, &
     217     &           ratqs,zqs,t)
     218
     219!===========================================================================
     220! Auteur : Arnaud Octavio Jam (LMD/CNRS)
     221! Date : 25 Mai 2010
     222! Objet : calcule les valeurs de qc et rneb dans les thermiques
     223!===========================================================================
     224
     225
     226      USE ioipsl_getin_p_mod, ONLY : getin_p
     227
     228      IMPLICIT NONE
     229
     230#include "YOMCST.h"
     231#include "YOETHF.h"
     232#include "FCTTRE.h"
     233#include "thermcell.h"
     234#include "nuage.h"
     235     
     236      INTEGER itap,ind1,ind2
     237      INTEGER ngrid,klev,klon,l,ig
     238     
     239      REAL ztv(ngrid,klev)
     240      REAL po(ngrid)
     241      REAL zqenv(ngrid)   
     242      REAL zqta(ngrid,klev)
     243         
     244      REAL fraca(ngrid,klev+1)
     245      REAL zpspsk(ngrid,klev)
     246      REAL paprs(ngrid,klev+1)
     247      REAL ztla(ngrid,klev)
     248      REAL zthl(ngrid,klev)
     249
     250      REAL zqsatth(ngrid,klev)
     251      REAL zqsatenv(ngrid,klev)
     252     
     253     
     254      REAL sigma1(ngrid,klev)                                                         
     255      REAL sigma2(ngrid,klev)
     256      REAL qlth(ngrid,klev)
     257      REAL qlenv(ngrid,klev)
     258      REAL qltot(ngrid,klev)
     259      REAL cth(ngrid,klev) 
     260      REAL cenv(ngrid,klev)   
     261      REAL ctot(ngrid,klev)
     262      REAL rneb(ngrid,klev)
     263      REAL t(ngrid,klev)                                                                 
     264      REAL qsatmmussig1,qsatmmussig2,sqrtpi,sqrt2,sqrt2pi,pi
     265      REAL rdd,cppd,Lv
     266      REAL alth,alenv,ath,aenv
     267      REAL sth,senv,sigma1s,sigma2s,xth,xenv,exp_xenv1,exp_xenv2,exp_xth1,exp_xth2
     268      REAL xth1,xth2,xenv1,xenv2,deltasth, deltasenv
     269      REAL IntJ,IntI1,IntI2,IntI3,coeffqlenv,coeffqlth
     270      REAL Tbef,zdelta,qsatbef,zcor
     271      REAL qlbef 
     272      REAL ratqs(ngrid,klev) ! determine la largeur de distribution de vapeur
     273      ! Change the width of the PDF used for vertical subgrid scale heterogeneity
     274      ! (J Jouhaud, JL Dufresne, JB Madeleine)
     275      REAL,SAVE :: vert_alpha
     276      LOGICAL, SAVE :: firstcall = .TRUE.
     277
     278      REAL zpdf_sig(ngrid),zpdf_k(ngrid),zpdf_delta(ngrid)
     279      REAL zpdf_a(ngrid),zpdf_b(ngrid),zpdf_e1(ngrid),zpdf_e2(ngrid)
     280      REAL zqs(ngrid), qcloud(ngrid)
     281      REAL erf
     282
     283!------------------------------------------------------------------------------
     284! Initialize
     285!------------------------------------------------------------------------------
     286      sigma1(:,:)=0.
     287      sigma2(:,:)=0.
     288      qlth(:,:)=0.
     289      qlenv(:,:)=0. 
     290      qltot(:,:)=0.
     291      rneb(:,:)=0.
     292      qcloud(:)=0.
     293      cth(:,:)=0.
     294      cenv(:,:)=0.
     295      ctot(:,:)=0.
     296      qsatmmussig1=0.
     297      qsatmmussig2=0.
     298      rdd=287.04
     299      cppd=1005.7
     300      pi=3.14159
     301      Lv=2.5e6
     302      sqrt2pi=sqrt(2.*pi)
     303      sqrt2=sqrt(2.)
     304      sqrtpi=sqrt(pi)
     305
     306      IF (firstcall) THEN
     307        vert_alpha=0.5
     308        CALL getin_p('cloudth_vert_alpha',vert_alpha)
     309        WRITE(*,*) 'cloudth_vert_alpha = ', vert_alpha
     310        firstcall=.FALSE.
     311      ENDIF
     312
     313!-------------------------------------------------------------------------------
     314! Calcul de la fraction du thermique et des ecart-types des distributions
     315!-------------------------------------------------------------------------------                 
     316      do ind1=1,ngrid
     317
     318      if ((ztv(ind1,1).gt.ztv(ind1,2)).and.(fraca(ind1,ind2).gt.1.e-10)) then !Thermal and environnement
     319
     320      zqenv(ind1)=(po(ind1)-fraca(ind1,ind2)*zqta(ind1,ind2))/(1.-fraca(ind1,ind2)) !qt = a*qtth + (1-a)*qtenv
     321
     322
     323      Tbef=zthl(ind1,ind2)*zpspsk(ind1,ind2)
     324      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
     325      qsatbef= R2ES*FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
     326      qsatbef=MIN(0.5,qsatbef)
     327      zcor=1./(1.-retv*qsatbef)
     328      qsatbef=qsatbef*zcor
     329      zqsatenv(ind1,ind2)=qsatbef
     330
     331
     332      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2)     !qsl, p84
     333      aenv=1./(1.+(alenv*Lv/cppd))                                      !al, p84
     334      senv=aenv*(po(ind1)-zqsatenv(ind1,ind2))                          !s, p84
     335
     336!zqenv = qt environnement
     337!zqsatenv = qsat environnement
     338!zthl = Tl environnement
     339
     340
     341      Tbef=ztla(ind1,ind2)*zpspsk(ind1,ind2)
     342      zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
     343      qsatbef= R2ES * FOEEW(Tbef,zdelta)/paprs(ind1,ind2)
     344      qsatbef=MIN(0.5,qsatbef)
     345      zcor=1./(1.-retv*qsatbef)
     346      qsatbef=qsatbef*zcor
     347      zqsatth(ind1,ind2)=qsatbef
     348           
     349      alth=(0.622*Lv*zqsatth(ind1,ind2))/(rdd*ztla(ind1,ind2)**2)       !qsl, p84
     350      ath=1./(1.+(alth*Lv/cppd))                                        !al, p84
     351      sth=ath*(zqta(ind1,ind2)-zqsatth(ind1,ind2))                      !s, p84
     352     
     353     
     354!zqta = qt thermals
     355!zqsatth = qsat thermals
     356!ztla = Tl thermals
     357
     358!------------------------------------------------------------------------------
     359! s standard deviation
     360!------------------------------------------------------------------------------
     361
     362      sigma1s=(1.1**0.5)*(fraca(ind1,ind2)**0.6)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
     363      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+0.002*zqta(ind1,ind2)
     364!      tests
     365!      sigma1s=(0.92**0.5)*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*((sth-senv)**2)**0.5+ratqs(ind1,ind2)*po(ind1)
     366!      sigma1s=(0.92*(fraca(ind1,ind2)**0.5)/(1-fraca(ind1,ind2))*(((sth-senv)**2)**0.5))+0.002*zqenv(ind1)
     367!      sigma2s=0.09*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.5+0.002*zqta(ind1,ind2)
     368!      sigma2s=(0.09*(((sth-senv)**2)**0.5)/((fraca(ind1,ind2)+0.02)**0.5))+ratqs(ind1,ind2)*zqta(ind1,ind2)
    392369!       if (paprs(ind1,ind2).gt.90000) then
    393370!       ratqs(ind1,ind2)=0.002
    … …  
    400377!      sigma2s=0.11*((sth-senv)**2)**0.5/(fraca(ind1,ind2)+0.02)**0.4+0.00003 
    401378 
    402 !------------------------------------------------------------------------------
    403 ! Calcul de l'eau condens?e et de la couverture nuageuse
    404 !------------------------------------------------------------------------------
    405       sqrt2pi=sqrt(2.*pi)
    406       xth=sth/(sqrt(2.)*sigma2s)
    407       xenv=senv/(sqrt(2.)*sigma1s)
    408       cth(ind1,ind2)=0.5*(1.+1.*erf(xth))
    409       cenv(ind1,ind2)=0.5*(1.+1.*erf(xenv))
    410       ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
    411 
    412       qlth(ind1,ind2)=sigma2s*((exp(-1.*xth**2)/sqrt2pi)+xth*sqrt(2.)*cth(ind1,ind2))
    413       qlenv(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))   
    414       qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
    415      
    416379       IF (iflag_cloudth_vert == 1) THEN
    417380!-------------------------------------------------------------------------------
    418 !  Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
    419 !-------------------------------------------------------------------------------
    420 !      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
    421 !      deltasth=ath*ratqs(ind1,ind2)*zqta(ind1,ind2)
     381!  Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
     382!-------------------------------------------------------------------------------
     383
    422384      deltasenv=aenv*ratqs(ind1,ind2)*zqsatenv(ind1,ind2)
    423385      deltasth=ath*ratqs(ind1,ind2)*zqsatth(ind1,ind2)
    424 !      deltasenv=aenv*0.01*po(ind1)
    425 !     deltasth=ath*0.01*zqta(ind1,ind2)   
     386
    426387      xenv1=(senv-deltasenv)/(sqrt(2.)*sigma1s)
    427388      xenv2=(senv+deltasenv)/(sqrt(2.)*sigma1s)
    … …  
    435396      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)   
    436397
     398      ! Environment
    437399      IntJ=sigma1s*(exp(-1.*xenv1**2)/sqrt2pi)+0.5*senv*(1+erf(xenv1))
    438400      IntI1=coeffqlenv*0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
    … …  
    441403
    442404      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
    443 !      qlenv(ind1,ind2)=IntJ
    444 !      print*, qlenv(ind1,ind2),'VERIF EAU'
    445 
    446 
     405
     406      ! Thermal
    447407      IntJ=sigma2s*(exp(-1.*xth1**2)/sqrt2pi)+0.5*sth*(1+erf(xth1))
    448 !      IntI1=coeffqlth*((0.5*xth1-xth2)*exp(-1.*xth1**2)+0.5*xth2*exp(-1.*xth2**2))
    449 !      IntI2=coeffqlth*0.5*sqrtpi*(0.5+xth2**2)*(erf(xth2)-erf(xth1))
    450408      IntI1=coeffqlth*0.5*(0.5*sqrtpi*(erf(xth2)-erf(xth1))+xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
    451409      IntI2=coeffqlth*xth2*(exp(-1.*xth2**2)-exp(-1.*xth1**2))
    452410      IntI3=coeffqlth*0.5*sqrtpi*xth2**2*(erf(xth2)-erf(xth1))
    453411      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
    454 !      qlth(ind1,ind2)=IntJ
    455 !      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
    456412      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
    457413
    … …  
    459415
    460416!-------------------------------------------------------------------------------
    461 !  Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
     417!  Version 3: Changes by J. Jouhaud; condensation for q > -delta s
    462418!-------------------------------------------------------------------------------
    463419!      deltasenv=aenv*ratqs(ind1,ind2)*po(ind1)
    … …  
    470426      xenv1=-(senv+deltasenv)/(sqrt(2.)*sigma1s)
    471427      xenv2=-(senv-deltasenv)/(sqrt(2.)*sigma1s)
     428      exp_xenv1 = exp(-1.*xenv1**2)
     429      exp_xenv2 = exp(-1.*xenv2**2)
    472430      xth1=-(sth+deltasth)/(sqrt(2.)*sigma2s)
    473431      xth2=-(sth-deltasth)/(sqrt(2.)*sigma2s)
    474 !     coeffqlenv=(sigma1s)**2/(2*sqrtpi*deltasenv)
    475 !     coeffqlth=(sigma2s)**2/(2*sqrtpi*deltasth)
     432      exp_xth1 = exp(-1.*xth1**2)
     433      exp_xth2 = exp(-1.*xth2**2)
    476434     
    477435      cth(ind1,ind2)=0.5*(1.-1.*erf(xth1))
    … …  
    479437      ctot(ind1,ind2)=fraca(ind1,ind2)*cth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*cenv(ind1,ind2)
    480438
    481       IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp(-1.*xenv2**2)
     439     
     440      !environnement
     441      IntJ=0.5*senv*(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
     442      if (deltasenv .lt. 1.e-10) then
     443      qlenv(ind1,ind2)=IntJ
     444      else
    482445      IntI1=(((senv+deltasenv)**2+(sigma1s)**2)/(8*deltasenv))*(erf(xenv2)-erf(xenv1))
    483       IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
    484       IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp(-1.*xenv1**2)-exp(-1.*xenv2**2))
    485 
    486 !      IntI1=0.5*(0.5*sqrtpi*(erf(xenv2)-erf(xenv1))+xenv1*exp(-1.*xenv1**2)-xenv2*exp(-1.*xenv2**2))
    487 !      IntI2=xenv2*(exp(-1.*xenv2**2)-exp(-1.*xenv1**2))
    488 !      IntI3=0.5*sqrtpi*xenv2**2*(erf(xenv2)-erf(xenv1))
    489 
     446      IntI2=(sigma1s**2/(4*deltasenv*sqrtpi))*(xenv1*exp_xenv1-xenv2*exp_xenv2)
     447      IntI3=((sqrt2*sigma1s*(senv+deltasenv))/(4*sqrtpi*deltasenv))*(exp_xenv1-exp_xenv2)
    490448      qlenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
    491 !      qlenv(ind1,ind2)=IntJ
    492 !      print*, qlenv(ind1,ind2),'VERIF EAU'
    493 
    494       IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp(-1.*xth2**2)
     449      endif
     450
     451
     452      !thermique
     453      IntJ=0.5*sth*(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
     454      if (deltasth .lt. 1.e-10) then ! Seuil a choisir !!!
     455      qlth(ind1,ind2)=IntJ
     456      else
    495457      IntI1=(((sth+deltasth)**2+(sigma2s)**2)/(8*deltasth))*(erf(xth2)-erf(xth1))
    496       IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp(-1.*xth1**2)-xth2*exp(-1.*xth2**2))
    497       IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp(-1.*xth1**2)-exp(-1.*xth2**2))
    498      
     458      IntI2=(sigma2s**2/(4*deltasth*sqrtpi))*(xth1*exp_xth1-xth2*exp_xth2)
     459      IntI3=((sqrt2*sigma2s*(sth+deltasth))/(4*sqrtpi*deltasth))*(exp_xth1-exp_xth2)
    499460      qlth(ind1,ind2)=IntJ+IntI1+IntI2+IntI3
    500 !      qlth(ind1,ind2)=IntJ
    501 !      print*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
     461      endif
     462
    502463      qltot(ind1,ind2)=fraca(ind1,ind2)*qlth(ind1,ind2)+(1.-1.*fraca(ind1,ind2))*qlenv(ind1,ind2)
    503      
    504 
    505464
    506465
    507466      ENDIF ! of if (iflag_cloudth_vert==1 or 2)
    508 
    509 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    510467
    511468      if (cenv(ind1,ind2).lt.1.e-10.or.cth(ind1,ind2).lt.1.e-10) then
    … …  
    515472      else
    516473               
    517       ctot(ind1,ind2)=ctot(ind1,ind2)
    518474      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqs(ind1)
    519475!      qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
    … …  
    521477
    522478      endif 
    523                        
    524      
    525          
    526 !     print*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
    527 
    528 
    529       else  ! gaussienne environnement seule
     479
     480      else  ! Environment only
    530481     
    531482      zqenv(ind1)=po(ind1)
    … …  
    539490     
    540491
    541 !      qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
     492!     qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
    542493      zthl(ind1,ind2)=t(ind1,ind2)*(101325/paprs(ind1,ind2))**(rdd/cppd)
    543494      alenv=(0.622*Lv*zqsatenv(ind1,ind2))/(rdd*zthl(ind1,ind2)**2) 
    … …  
    548499      sigma1s=ratqs(ind1,ind2)*zqenv(ind1)
    549500
    550       sqrt2pi=sqrt(2.*pi)
    551       xenv=senv/(sqrt(2.)*sigma1s)
     501      xenv=senv/(sqrt2*sigma1s)
    552502      ctot(ind1,ind2)=0.5*(1.+1.*erf(xenv))
    553       qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt(2.)*cenv(ind1,ind2))
     503      qltot(ind1,ind2)=sigma1s*((exp(-1.*xenv**2)/sqrt2pi)+xenv*sqrt2*cenv(ind1,ind2))
    554504     
    555505      if (ctot(ind1,ind2).lt.1.e-3) then
    … …  
    559509      else   
    560510               
    561       ctot(ind1,ind2)=ctot(ind1,ind2)
    562511      qcloud(ind1)=qltot(ind1,ind2)/ctot(ind1,ind2)+zqsatenv(ind1,ind2)
    563512
    564513      endif   
    565514 
    566  
    567  
    568  
    569  
    570  
    571       endif   
    572       enddo
     515      endif       ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
     516      enddo       ! from the loop on ngrid l.333
    573517     
    574518      return

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

    r2594 r2669  
    1919       ok_ade, ok_aie, ok_cdnc, aerosol_couple, &
    2020       flag_aerosol, flag_aerosol_strat, new_aod, &
    21        bl95_b0, bl95_b1,&
     21       flag_bc_internal_mixture, bl95_b0, bl95_b1,&
    2222       read_climoz, &
    2323       alp_offset)
    … …  
    6464    ! ok_cdnc, ok cloud droplet number concentration
    6565    ! flag_aerosol_strat : flag pour les aerosols stratos
     66    ! flag_bc_internal_mixture : use BC internal mixture if true
    6667    ! bl95_b*: parameters in the formula to link CDNC to aerosol mass conc
    6768    !
    … …  
    7778    INTEGER              :: flag_aerosol
    7879    INTEGER              :: flag_aerosol_strat
     80    LOGICAL              :: flag_bc_internal_mixture
    7981    LOGICAL              :: new_aod
    8082    REAL                 :: bl95_b0, bl95_b1
    … …  
    9597    INTEGER, SAVE       :: flag_aerosol_omp
    9698    INTEGER, SAVE       :: flag_aerosol_strat_omp
     99    LOGICAL, SAVE       :: flag_bc_internal_mixture_omp
    97100    LOGICAL, SAVE       :: new_aod_omp
    98101    REAL,SAVE           :: bl95_b0_omp, bl95_b1_omp
    … …  
    398401    flag_aerosol_strat_omp = 0
    399402    CALL getin('flag_aerosol_strat',flag_aerosol_strat_omp)
     403
     404    !
     405    !Config Key  = flag_bc_internal_mixture
     406    !Config Desc = state of mixture for BC aerosols
     407    ! - n = external mixture
     408    ! - y = internal mixture
     409    !Config Def  = n
     410    !Config Help = Used in physiq.F / aeropt
     411    !
     412    flag_bc_internal_mixture_omp = .false.
     413    CALL getin('flag_bc_internal_mixture',flag_bc_internal_mixture_omp)
    400414
    401415    ! Temporary variable for testing purpose!
    … …  
    21362150    flag_aerosol=flag_aerosol_omp
    21372151    flag_aerosol_strat=flag_aerosol_strat_omp
     2152    flag_bc_internal_mixture=flag_bc_internal_mixture_omp
    21382153    new_aod=new_aod_omp
    21392154    aer_type = aer_type_omp
    … …  
    23112326    IF (ok_aie .AND. .NOT. ok_cdnc) THEN
    23122327       CALL abort_physic('conf_phys', 'ok_cdnc must be set to y if ok_aie is activated',1)
     2328    ENDIF
     2329
     2330    ! BC internal mixture is only possible with RRTM & NSW=6 & flag_aerosol=6 or aerosol_couple
     2331    IF (flag_bc_internal_mixture .AND. NSW.NE.6) THEN
     2332       CALL abort_physic('conf_phys', 'flag_bc_internal_mixture can only be activated with NSW=6',1)
     2333    ENDIF
     2334    IF (flag_bc_internal_mixture .AND. iflag_rrtm.NE.1) THEN
     2335       CALL abort_physic('conf_phys', 'flag_bc_internal_mixture can only be activated with RRTM',1)
     2336    ENDIF
     2337    IF (flag_bc_internal_mixture .AND. flag_aerosol.NE.6) THEN
     2338       CALL abort_physic('conf_phys', 'flag_bc_internal_mixture can only be activated with flag_aerosol=6',1)
    23132339    ENDIF
    23142340

  • LMDZ5/branches/testing/libf/phylmd/cosp/cosp_output_mod.F90

    r2594 r2669  
    4242         "clmcalipso", "Lidar Mid-level Cloud Fraction", "%", (/ ('', i=1, 3) /))
    4343  TYPE(ctrl_outcosp), SAVE :: o_clhcalipso = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
    44          "clhcalipso", "Lidar Hight-level Cloud Fraction", "%", (/ ('', i=1, 3) /))
     44         "clhcalipso", "Lidar High-level Cloud Fraction", "%", (/ ('', i=1, 3) /))
    4545  TYPE(ctrl_outcosp), SAVE :: o_cltcalipso = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
    4646         "cltcalipso", "Lidar Total Cloud Fraction", "%", (/ ('', i=1, 3) /))
    … …  
    145145         "clmmodis", "MODIS Mid-level Cloud Fraction", "1", (/ ('', i=1, 3) /))
    146146  TYPE(ctrl_outcosp), SAVE :: o_clhmodis = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
    147          "clhmodis", "MODIS Hight-level Cloud Fraction", "1", (/ ('', i=1, 3) /))
     147         "clhmodis", "MODIS High-level Cloud Fraction", "1", (/ ('', i=1, 3) /))
    148148  TYPE(ctrl_outcosp), SAVE :: o_cltmodis = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
    149149         "cltmodis", "MODIS Total Cloud Fraction", "1", (/ ('', i=1, 3) /))

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

    r2641 r2669  
    10771077        DO k = 1,nd
    10781078        write (6, '(i4,5(1x,e13.6))'), &
    1079           k, mp(idcum(igout),k), water(idcum(igout),k), ice(idcum(igout),k), &
    1080            evap(idcum(igout),k), fondue(idcum(igout),k)
     1079          k, mp(igout,k), water(igout,k), ice(igout,k), &
     1080           evap(igout,k), fondue(igout,k)
    10811081        ENDDO
    10821082        Print *, 'cva_driver after cv3_unsat: wdtrainA, wdtrainM '
    10831083        DO k = 1,nd
    10841084        write (6, '(i4,2(1x,e13.6))'), &
    1085            k, wdtrainA(idcum(igout),k), wdtrainM(idcum(igout),k)
     1085           k, wdtrainA(igout,k), wdtrainM(igout,k)
    10861086        ENDDO
    10871087      ENDIF
    … …  
    11291129!
    11301130      IF (debut) THEN
    1131         PRINT *, ' cv3_yield -> fqd(1) = ', fqd(idcum(igout), 1)
     1131        PRINT *, ' cv3_yield -> fqd(1) = ', fqd(igout, 1)
    11321132      END IF !(debut) THEN
    11331133!   
    11341134      IF (prt_level >= 10) THEN
    11351135        Print *, 'cva_driver after cv3_yield:ft(1) , ftd(1) ', &
    1136                     ft(idcum(igout),1), ftd(idcum(igout),1)
     1136                    ft(igout,1), ftd(igout,1)
    11371137        Print *, 'cva_driver after cv3_yield:fq(1) , fqd(1) ', &
    1138                     fq(idcum(igout),1), fqd(idcum(igout),1)
     1138                    fq(igout,1), fqd(igout,1)
    11391139      ENDIF
    11401140!   

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

    r2073 r2669  
    99  INTEGER,SAVE :: klevm1
    1010  INTEGER,SAVE :: kflev
    11   INTEGER,SAVE :: nslay
    1211
    13 !$OMP THREADPRIVATE(klon,kfdia,kidia,kdlon,nslay)
     12!$OMP THREADPRIVATE(klon,kfdia,kidia,kdlon)
    1413  REAL,save,allocatable,dimension(:) :: zmasq
    1514!$OMP THREADPRIVATE(zmasq)   
    … …  
    2423   
    2524    klon=klon0
    26     nslay=1 ! Slab, provisoire (F. Codron)
    2725    kdlon=klon
    2826    kidia=1

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

    r2542 r2669  
    624624                    zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i)))
    625625                    else if (iflag_t_glace.ge.1) then
    626                     zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
     626                    zdelta = MAX(0.,SIGN(1.,t_glace_max-Tbef(i)))
     627! BUG corrige par JYG   zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i)))
    627628                    endif
    628629                 endif

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

    r2542 r2669  
    1717#ifdef CPP_XIOS
    1818  INTERFACE histwrite_phy
    19     MODULE PROCEDURE histwrite2d_phy,histwrite3d_phy,histwrite2d_phy_old,histwrite3d_phy_old,histwrite2d_xios,histwrite3d_xios
     19!#ifdef CPP_XIOSnew
     20    MODULE PROCEDURE histwrite2d_phy,histwrite3d_phy,histwrite2d_phy_old,histwrite3d_phy_old,histwrite2d_xios,histwrite3d_xios,histwrite0d_xios
     21!#else
     22!    MODULE PROCEDURE histwrite2d_phy,histwrite3d_phy,histwrite2d_phy_old,histwrite3d_phy_old,histwrite2d_xios,histwrite3d_xios
     23!#endif
     24
    2025  END INTERFACE
    2126#else
    … …  
    13311336  IF (prt_level >= 10) write(lunout,*)'End histrwrite3d_xios ',field_name
    13321337  END SUBROUTINE histwrite3d_xios
     1338
     1339#ifdef CPP_XIOS
     1340  SUBROUTINE histwrite0d_xios(field_name, field)
     1341  USE xios, only: xios_send_field
     1342  IMPLICIT NONE
     1343
     1344    CHARACTER(LEN=*), INTENT(IN) :: field_name
     1345    REAL, INTENT(IN) :: field ! --> scalar
     1346
     1347!$OMP MASTER
     1348   CALL xios_send_field(field_name, field)
     1349!$OMP END MASTER
     1350
     1351  END SUBROUTINE histwrite0d_xios
     1352#endif
     1353
    13331354#endif
    13341355end module iophy

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

    r2542 r2669  
    164164  diff_sst(:) = diff_sst_save(:)
    165165  diff_siv(:) = diff_siv_save(:)
    166 ! For Debug purpose
    167 !  PRINT *,'limit_slab sst',MINVAL(diff_sst(:)),MAXVAL(diff_sst(:))
    168 !  PRINT *,'limit_slab siv',MINVAL(diff_siv(:)),MAXVAL(diff_siv(:))
    169 !  PRINT *,'limit_slab bils',MINVAL(lmt_bils(:)),MAXVAL(lmt_bils(:))
    170166
    171167END SUBROUTINE limit_slab

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

    r2408 r2669  
    99  USE indice_sol_mod
    1010  USE surface_data
     11  USE mod_grid_phy_lmdz, ONLY: klon_glo
     12  USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
    1113
    1214  IMPLICIT NONE
    … …  
    1416  PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice, ocean_slab_ice
    1517
    16 !****************************************************************************************
     18!***********************************************************************************
    1719! Global saved variables
    18 !****************************************************************************************
    19 
     20!***********************************************************************************
     21  ! number of slab vertical layers
     22  INTEGER, PUBLIC, SAVE :: nslay
     23  !$OMP THREADPRIVATE(nslay)
     24  ! timestep for coupling (update slab temperature) in timesteps
    2025  INTEGER, PRIVATE, SAVE                           :: cpl_pas
    2126  !$OMP THREADPRIVATE(cpl_pas)
     27  ! cyang = 1/heat capacity of top layer (rho.c.H)
    2228  REAL, PRIVATE, SAVE                              :: cyang
    2329  !$OMP THREADPRIVATE(cyang)
     30  ! depth of slab layers (1 or 2)
    2431  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: slabh
    2532  !$OMP THREADPRIVATE(slabh)
     33  ! slab temperature
    2634  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: tslab
    2735  !$OMP THREADPRIVATE(tslab)
     36  ! heat flux convergence due to Ekman
     37  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_ekman
     38  !$OMP THREADPRIVATE(dt_ekman)
     39  ! heat flux convergence due to horiz diffusion
     40  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_hdiff
     41  !$OMP THREADPRIVATE(dt_hdiff)
     42  ! fraction of ocean covered by sea ice (sic / (oce+sic))
    2843  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: fsic
    2944  !$OMP THREADPRIVATE(fsic)
     45  ! temperature of the sea ice
    3046  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: tice
    3147  !$OMP THREADPRIVATE(tice)
     48  ! sea ice thickness, in kg/m2
    3249  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: seaice
    3350  !$OMP THREADPRIVATE(seaice)
     51  ! net surface heat flux, weighted by open ocean fraction
    3452  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bils
    3553  !$OMP THREADPRIVATE(slab_bils)
     54  ! slab_bils accumulated over cpl_pas timesteps
    3655  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bils_cum
    3756  !$OMP THREADPRIVATE(bils_cum)
     57  ! net heat flux into the ocean below the ice : conduction + solar radiation
    3858  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bilg
    3959  !$OMP THREADPRIVATE(slab_bilg)
     60  ! slab_bilg over cpl_pas timesteps
    4061  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bilg_cum
    4162  !$OMP THREADPRIVATE(bilg_cum)
    42 
    43 !****************************************************************************************
     63  ! wind stress saved over cpl_pas timesteps
     64  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: taux_cum
     65  !$OMP THREADPRIVATE(taux_cum)
     66  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: tauy_cum
     67  !$OMP THREADPRIVATE(tauy_cum)
     68
     69!***********************************************************************************
    4470! Parameters (could be read in def file: move to slab_init)
    45 !****************************************************************************************
     71!***********************************************************************************
    4672! snow and ice physical characteristics:
    4773    REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp
    4874    REAL, PARAMETER :: t_melt=273.15   ! melting ice temp
    4975    REAL, PARAMETER :: sno_den=300. !mean snow density, kg/m3
    50     REAL, PARAMETER :: ice_den=917.
    51     REAL, PARAMETER :: sea_den=1025.
    52     REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity
    53     REAL, PARAMETER :: sno_cond=0.31*sno_den
     76    REAL, PARAMETER :: ice_den=917. ! ice density
     77    REAL, PARAMETER :: sea_den=1025. ! sea water density
     78    REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity of ice
     79    REAL, PARAMETER :: sno_cond=0.31*sno_den ! conductivity of snow
    5480    REAL, PARAMETER :: ice_cap=2067.   ! specific heat capacity, snow and ice
     81    REAL, PARAMETER :: sea_cap=3995.   ! specific heat capacity, snow and ice
    5582    REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice
    5683
    … …  
    74101    REAL, PARAMETER :: alb_ice_dry=0.75 !dry thick ice
    75102    REAL, PARAMETER :: alb_ice_wet=0.66 !melting thick ice
    76     REAL, PARAMETER :: pen_frac=0.3 !fraction of penetrating shortwave (no snow)
     103    REAL, PARAMETER :: pen_frac=0.3 !fraction of shortwave penetrating into the
     104    ! ice (no snow)
    77105    REAL, PARAMETER :: pen_ext=1.5 !extinction of penetrating shortwave (m-1)
    78106
    79 !****************************************************************************************
     107! horizontal transport
     108   LOGICAL, PUBLIC, SAVE :: slab_hdiff
     109   !$OMP THREADPRIVATE(slab_hdiff)
     110   REAL, PRIVATE, SAVE    :: coef_hdiff ! coefficient for horizontal diffusion
     111   !$OMP THREADPRIVATE(coef_hdiff)
     112   INTEGER, PUBLIC, SAVE :: slab_ekman, slab_cadj ! Ekman, conv adjustment
     113   !$OMP THREADPRIVATE(slab_ekman)
     114   !$OMP THREADPRIVATE(slab_cadj)
     115
     116!***********************************************************************************
    80117
    81118CONTAINS
    82119!
    83 !****************************************************************************************
     120!***********************************************************************************
    84121!
    85122  SUBROUTINE ocean_slab_init(dtime, pctsrf_rst)
    86123  !, seaice_rst etc
    87124
    88     use IOIPSL
    89 
    90     ! For ok_xxx vars (Ekman...)
    91     INCLUDE "clesphys.h"
     125    USE ioipsl_getin_p_mod, ONLY : getin_p
     126    USE mod_phys_lmdz_transfert_para, ONLY : gather
     127    USE slab_heat_transp_mod, ONLY : ini_slab_transp
    92128
    93129    ! Input variables
    94 !****************************************************************************************
     130!***********************************************************************************
    95131    REAL, INTENT(IN)                         :: dtime
    96132! Variables read from restart file
    97     REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
     133    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
     134    ! surface fractions from start file
    98135
    99136! Local variables
    100 !****************************************************************************************
     137!************************************************************************************
    101138    INTEGER                :: error
     139    REAL, DIMENSION(klon_glo) :: zmasq_glo
    102140    CHARACTER (len = 80)   :: abort_message
    103141    CHARACTER (len = 20)   :: modname = 'ocean_slab_intit'
    104142
    105 !****************************************************************************************
     143!***********************************************************************************
     144! Define some parameters
     145!***********************************************************************************
     146! Number of slab layers
     147    nslay=2
     148    CALL getin_p('slab_layers',nslay)
     149    print *,'number of slab layers : ',nslay
     150! Layer thickness
     151    ALLOCATE(slabh(nslay), stat = error)
     152    IF (error /= 0) THEN
     153       abort_message='Pb allocation slabh'
     154       CALL abort_physic(modname,abort_message,1)
     155    ENDIF
     156    slabh(1)=50.
     157    IF (nslay.GT.1) THEN
     158        slabh(2)=150.
     159    END IF
     160
     161! cyang = 1/heat capacity of top layer (rho.c.H)
     162    cyang=1/(slabh(1)*sea_den*sea_cap)
     163
     164! cpl_pas  coupling period (update of tslab and ice fraction)
     165! pour un calcul a chaque pas de temps, cpl_pas=1
     166    cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
     167    CALL getin_p('cpl_pas',cpl_pas)
     168    print *,'cpl_pas',cpl_pas
     169
     170! Horizontal diffusion
     171    slab_hdiff=.FALSE.
     172    CALL getin_p('slab_hdiff',slab_hdiff)
     173    coef_hdiff=25000.
     174    CALL getin_p('coef_hdiff',coef_hdiff)
     175! Ekman transport
     176    slab_ekman=0
     177    CALL getin_p('slab_ekman',slab_ekman)
     178! Convective adjustment
     179    IF (nslay.EQ.1) THEN
     180        slab_cadj=0
     181    ELSE
     182        slab_cadj=1
     183    END IF
     184    CALL getin_p('slab_cadj',slab_cadj)
     185
     186!************************************************************************************
    106187! Allocate surface fraction read from restart file
    107 !****************************************************************************************
     188!************************************************************************************
    108189    ALLOCATE(fsic(klon), stat = error)
    109190    IF (error /= 0) THEN
    … …  
    112193    ENDIF
    113194    fsic(:)=0.
     195    !zmasq = continent fraction
    114196    WHERE (1.-zmasq(:)>EPSFRA)
    115197        fsic(:) = pctsrf_rst(:,is_sic)/(1.-zmasq(:))
    116198    END WHERE
    117199
    118 !****************************************************************************************
    119 ! Allocate saved variables
    120 !****************************************************************************************
     200!************************************************************************************
     201! Allocate saved fields
     202!************************************************************************************
    121203    ALLOCATE(tslab(klon,nslay), stat=error)
    122204       IF (error /= 0) CALL abort_physic &
    … …  
    136218    bils_cum(:) = 0.0   
    137219
    138     IF (version_ocean=='sicINT') THEN
     220    IF (version_ocean=='sicINT') THEN ! interactive sea ice
    139221        ALLOCATE(slab_bilg(klon), stat = error)
    140222        IF (error /= 0) THEN
    … …  
    161243    END IF
    162244
    163 !****************************************************************************************
    164 ! Define some parameters
    165 !****************************************************************************************
    166 ! Layer thickness
    167     ALLOCATE(slabh(nslay), stat = error)
    168     IF (error /= 0) THEN
    169        abort_message='Pb allocation slabh'
    170        CALL abort_physic(modname,abort_message,1)
     245    IF (slab_hdiff) THEN !horizontal diffusion
     246        ALLOCATE(dt_hdiff(klon,nslay), stat = error)
     247        IF (error /= 0) THEN
     248           abort_message='Pb allocation dt_hdiff'
     249           CALL abort_physic(modname,abort_message,1)
     250        ENDIF
     251        dt_hdiff(:,:) = 0.0   
    171252    ENDIF
    172     slabh(1)=50.
    173 ! cyang = 1/heat capacity of top layer (rho.c.H)
    174     cyang=1/(slabh(1)*4.228e+06)
    175 
    176 ! cpl_pas periode de couplage avec slab (update tslab, pctsrf)
    177 ! pour un calcul à chaque pas de temps, cpl_pas=1
    178     cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
    179     CALL getin('cpl_pas',cpl_pas)
    180     print *,'cpl_pas',cpl_pas
     253
     254    IF (slab_ekman.GT.0) THEN ! ekman transport
     255        ALLOCATE(dt_ekman(klon,nslay), stat = error)
     256        IF (error /= 0) THEN
     257           abort_message='Pb allocation dt_ekman'
     258           CALL abort_physic(modname,abort_message,1)
     259        ENDIF
     260        dt_ekman(:,:) = 0.0   
     261        ALLOCATE(taux_cum(klon), stat = error)
     262        IF (error /= 0) THEN
     263           abort_message='Pb allocation taux_cum'
     264           CALL abort_physic(modname,abort_message,1)
     265        ENDIF
     266        taux_cum(:) = 0.0   
     267        ALLOCATE(tauy_cum(klon), stat = error)
     268        IF (error /= 0) THEN
     269           abort_message='Pb allocation tauy_cum'
     270           CALL abort_physic(modname,abort_message,1)
     271        ENDIF
     272        tauy_cum(:) = 0.0   
     273    ENDIF
     274
     275! Initialize transport
     276    IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN
     277      CALL gather(zmasq,zmasq_glo)
     278!$OMP MASTER  ! Only master thread
     279      IF (is_mpi_root) THEN
     280          CALL ini_slab_transp(zmasq_glo)
     281      END IF
     282!$OMP END MASTER
     283    END IF
    181284
    182285 END SUBROUTINE ocean_slab_init
    183286!
    184 !****************************************************************************************
     287!***********************************************************************************
    185288!
    186289  SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified)
    187290
    188     USE limit_read_mod
    189 
    190 !    INCLUDE "clesphys.h"
     291! this routine sends back the sea ice and ocean fraction to the main physics
     292! routine. Called only with interactive sea ice
    191293
    192294! Arguments
    193 !****************************************************************************************
    194     INTEGER, INTENT(IN)                        :: itime   ! numero du pas de temps courant
    195     INTEGER, INTENT(IN)                        :: jour    ! jour a lire dans l'annee
    196     REAL   , INTENT(IN)                        :: dtime   ! pas de temps de la physique (en s)
     295!************************************************************************************
     296    INTEGER, INTENT(IN)                        :: itime   ! current timestep
     297    INTEGER, INTENT(IN)                        :: jour    !  day in year (not
     298    REAL   , INTENT(IN)                        :: dtime   ! physics timestep (s)
    197299    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg  ! sub-surface fraction
    198     LOGICAL, INTENT(OUT)                       :: is_modified ! true if pctsrf is modified at this time step
    199 
    200 ! Local variables
    201 !****************************************************************************************
    202 
    203     IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN   
    204        CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, is_modified)
    205     ELSE
     300    LOGICAL, INTENT(OUT)                       :: is_modified ! true if pctsrf is
     301                                                         ! modified at this time step
     302
    206303       pctsrf_chg(:,is_oce)=(1.-fsic(:))*(1.-zmasq(:))
    207304       pctsrf_chg(:,is_sic)=fsic(:)*(1.-zmasq(:))
    208305       is_modified=.TRUE.
    209     END IF
    210306
    211307  END SUBROUTINE ocean_slab_frac
    212308!
    213 !****************************************************************************************
     309!************************************************************************************
    214310!
    215311  SUBROUTINE ocean_slab_noice( &
    … …  
    224320   
    225321    USE calcul_fluxs_mod
     322    USE slab_heat_transp_mod, ONLY: divgrad_phy, slab_ekman1, slab_ekman2
     323    USE mod_phys_lmdz_para
    226324
    227325    INCLUDE "clesphys.h"
    228326
     327! This routine
     328! (1) computes surface turbulent fluxes over points with some open ocean   
     329! (2) reads additional Q-flux (everywhere)
     330! (3) computes horizontal transport (diffusion & Ekman)
     331! (4) updates slab temperature every cpl_pas ; creates new ice if needed.
     332
     333! Note :
     334! klon total number of points
     335! knon number of points with open ocean (varies with time)
     336! knindex gives position of the knon points within klon.
     337! In general, local saved variables have klon values
     338! variables exchanged with PBL module have knon.
     339
    229340! Input arguments
    230 !****************************************************************************************
    231     INTEGER, INTENT(IN)                  :: itime
    232     INTEGER, INTENT(IN)                  :: jour
    233     INTEGER, INTENT(IN)                  :: knon
    234     INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
    235     REAL, INTENT(IN)                     :: dtime
    236     REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
    237     REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragq, cdragm
    238     REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
    239     REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
    240     REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
    241     REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
    242     REAL, DIMENSION(klon), INTENT(IN)    :: ps
    243     REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness
    244     REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
    245     REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
     341!***********************************************************************************
     342    INTEGER, INTENT(IN)                  :: itime ! current timestep INTEGER,
     343    INTEGER, INTENT(IN)                  :: jour  ! day in year (for Q-Flux)
     344    INTEGER, INTENT(IN)                  :: knon  ! number of points
     345    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
     346    REAL, INTENT(IN) :: dtime  ! timestep (s)
     347    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
     348    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragq, cdragm
     349    ! drag coefficients
     350    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
     351    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum ! near surface T, q
     352    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
     353    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
     354    ! exchange coefficients for boundary layer scheme
     355    REAL, DIMENSION(klon), INTENT(IN)    :: ps  ! surface pressure
     356    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness ! surface wind
     357    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in ! surface temperature
     358    REAL, DIMENSION(klon), INTENT(INOUT) :: radsol ! net surface radiative flux
    246359
    247360! In/Output arguments
    248 !****************************************************************************************
    249     REAL, DIMENSION(klon), INTENT(INOUT) :: snow
     361!************************************************************************************
     362    REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! in kg/m2
    250363   
    251364! Output arguments
    252 !****************************************************************************************
     365!************************************************************************************
    253366    REAL, DIMENSION(klon), INTENT(OUT)   :: qsurf
    254367    REAL, DIMENSION(klon), INTENT(OUT)   :: evap, fluxsens, fluxlat
    255368    REAL, DIMENSION(klon), INTENT(OUT)   :: flux_u1, flux_v1
    256     REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new
     369    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new ! new surface tempearture
    257370    REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l     
    258371    REAL, DIMENSION(klon), INTENT(OUT)   :: qflux
    259372
    260373! Local variables
    261 !****************************************************************************************
    262     INTEGER               :: i,ki
    263     ! surface fluxes
     374!************************************************************************************
     375    INTEGER               :: i,ki,k
     376    REAL                  :: t_cadj
     377    !  for surface heat fluxes
    264378    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
    265     ! flux correction
     379    ! for Q-Flux computation: d/dt SST, d/dt ice volume (kg/m2), surf fluxes
    266380    REAL, DIMENSION(klon) :: diff_sst, diff_siv, lmt_bils
    267     ! surface wind stress
     381    ! for surface wind stress
    268382    REAL, DIMENSION(klon) :: u0, v0
    269383    REAL, DIMENSION(klon) :: u1_lay, v1_lay
    270     ! ice creation
     384    ! for new ice creation
    271385    REAL                  :: e_freeze, h_new, dfsic
    272 
    273 !****************************************************************************************
    274 ! 1) Flux calculation
    275 !
    276 !****************************************************************************************
    277     cal(:)      = 0.
    278     beta(:)     = 1.
    279     dif_grnd(:) = 0.
     386    ! horizontal diffusion and Ekman local vars
     387    ! dimension = global domain (klon_glo) instead of // subdomains
     388    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo, dt_ekman_glo
     389    ! dt_ekman_glo saved for diagnostic, dt_ekman_tmp used for time loop
     390    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_tmp, dt_ekman_tmp
     391    REAL, DIMENSION(klon_glo,nslay) :: tslab_glo
     392    REAL, DIMENSION(klon_glo) :: taux_glo,tauy_glo
     393
     394!****************************************************************************************
     395! 1) Surface fluxes calculation
     396!   
     397!****************************************************************************************
     398    cal(:)      = 0. ! infinite thermal inertia
     399    beta(:)     = 1. ! wet surface
     400    dif_grnd(:) = 0. ! no diffusion into ground
    280401   
    281402! Suppose zero surface speed
    … …  
    285406    v1_lay(:) = v1(:) - v0(:)
    286407
     408! Compute latent & sensible fluxes
    287409    CALL calcul_fluxs(knon, is_oce, dtime, &
    288410         tsurf_in, p1lay, cal, beta, cdragh, cdragq, ps, &
    … …  
    292414         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
    293415
    294 ! - Flux calculation at first modele level for U and V
    295     CALL calcul_flux_wind(knon, dtime, &
    296          u0, v0, u1, v1, gustiness, cdragm, &
    297          AcoefU, AcoefV, BcoefU, BcoefV, &
    298          p1lay, temp_air, &
    299          flux_u1, flux_v1) 
    300 
    301 ! Accumulate total fluxes locally
     416! save total cumulated heat fluxes locally
     417! radiative + turbulent + melt of falling snow
    302418    slab_bils(:)=0.
    303419    DO i=1,knon
    … …  
    306422                      -precip_snow(i)*ice_lat*(1.+snow_wfact*fsic(ki)))
    307423        bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
    308 ! Also taux, tauy, saved vars...
    309424    END DO
    310425
    311 !****************************************************************************************
    312 ! 2) Get global variables lmt_bils and diff_sst from file limit_slab.nc
     426!  Compute surface wind stress
     427    CALL calcul_flux_wind(knon, dtime, &
     428         u0, v0, u1, v1, gustiness, cdragm, &
     429         AcoefU, AcoefV, BcoefU, BcoefV, &
     430         p1lay, temp_air, &
     431         flux_u1, flux_v1) 
     432
     433! save cumulated wind stress
     434    IF (slab_ekman.GT.0) THEN
     435      DO i=1,knon
     436          ki=knindex(i)
     437          taux_cum(ki)=taux_cum(ki)+flux_u1(i)*(1.-fsic(ki))/cpl_pas
     438          tauy_cum(ki)=tauy_cum(ki)+flux_v1(i)*(1.-fsic(ki))/cpl_pas
     439      END DO
     440    ENDIF
     441
     442!****************************************************************************************
     443! 2) Q-Flux : get global variables lmt_bils, diff_sst and diff_siv from file limit_slab.nc
    313444!
    314445!****************************************************************************************
    315446    lmt_bils(:)=0.
    316     CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv) ! global pour un processus
     447    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv)
    317448    ! lmt_bils and diff_sst,siv saved by limit_slab
    318449    qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.-diff_siv(:)*ice_den*ice_lat/86400.
    … …  
    320451
    321452!****************************************************************************************
    322 ! 3) Recalculate new temperature
    323 !
     453! 3) Recalculate new temperature (add Surf fluxes, Q-Flux, Ocean transport)
     454!    Bring to freezing temp and make sea ice if necessary
     455! 
    324456!***********************************************o*****************************************
    325457    tsurf_new=tsurf_in
    326458    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
    327         ! Compute transport
    328         ! Add read QFlux and SST tendency
    329         tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas
    330         ! Add cumulated surface fluxes
    331         tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
    332         ! Update surface temperature
    333         SELECT CASE(version_ocean)
    334         CASE('sicNO')
    335             DO i=1,knon
    336                 ki=knindex(i)
    337                 tsurf_new(i)=tslab(ki,1)
     459! ***********************************
     460!  Horizontal transport
     461! ***********************************
     462      IF (slab_ekman.GT.0) THEN
     463          ! copy wind stress to global var
     464          CALL gather(taux_cum,taux_glo)
     465          CALL gather(tauy_cum,tauy_glo)
     466      END IF
     467
     468      IF (slab_hdiff.OR.(slab_ekman.GT.0)) THEN
     469        CALL gather(tslab,tslab_glo)
     470      ! Compute horiz transport on one process only
     471!$OMP MASTER  ! Only master thread
     472        IF (is_mpi_root) THEN ! Only master processus         
     473          IF (slab_hdiff) THEN
     474              dt_hdiff_glo(:,:)=0.
     475          END IF
     476          IF (slab_ekman.GT.0) THEN
     477              dt_ekman_glo(:,:)=0.
     478          END IF
     479          DO i=1,cpl_pas ! time splitting for numerical stability
     480            IF (slab_ekman.GT.0) THEN
     481              SELECT CASE (slab_ekman)
     482                CASE (1)
     483                  CALL slab_ekman1(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
     484                CASE (2)
     485                  CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
     486                CASE DEFAULT
     487                  dt_ekman_tmp(:,:)=0.
     488              END SELECT
     489              dt_ekman_glo(:,:)=dt_ekman_glo(:,:)+dt_ekman_tmp(:,:)
     490              ! convert dt_ekman from K.s-1.(kg.m-2) to K.s-1   
     491              DO k=1,nslay
     492                dt_ekman_tmp(:,k)=dt_ekman_tmp(:,k)/(slabh(k)*sea_den)
     493              ENDDO
     494              tslab_glo=tslab_glo+dt_ekman_tmp*dtime
     495            ENDIF
     496            IF (slab_hdiff) THEN ! horizontal diffusion
     497              ! laplacian of slab T
     498              CALL divgrad_phy(nslay,tslab_glo,dt_hdiff_tmp)
     499              ! multiply by diff coef and normalize to 50m slab equivalent
     500              dt_hdiff_tmp=dt_hdiff_tmp*coef_hdiff*50./SUM(slabh)
     501              dt_hdiff_glo(:,:)=dt_hdiff_glo(:,:)+ dt_hdiff_tmp(:,:)
     502              tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
     503            END IF
     504          END DO ! time splitting
     505          IF (slab_hdiff) THEN
     506            !dt_hdiff_glo saved in W/m2
     507            DO k=1,nslay
     508              dt_hdiff_glo(:,k)=dt_hdiff_glo(:,k)*slabh(k)*sea_den*sea_cap/cpl_pas
    338509            END DO
    339         CASE('sicOBS') ! check for sea ice or tslab below freezing
    340             DO i=1,knon
    341                 ki=knindex(i)
    342                 IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
    343                     tslab(ki,1)=t_freeze
    344                 END IF
    345                 tsurf_new(i)=tslab(ki,1)
    346             END DO
    347         CASE('sicINT')
    348             DO i=1,knon
    349                 ki=knindex(i)
    350                 IF (fsic(ki).LT.epsfra) THEN ! Free of ice
    351                     IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
    352                         ! quantity of new ice formed
    353                         e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
    354                         ! new ice
    355                         tice(ki)=t_freeze
    356                         fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
    357                         IF (fsic(ki).GT.ice_frac_min) THEN
    358                             seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
    359                             tslab(ki,1)=t_freeze
    360                         ELSE
    361                             fsic(ki)=0.
    362                         END IF
    363                         tsurf_new(i)=t_freeze
    364                     ELSE
    365                         tsurf_new(i)=tslab(ki,1)
    366                     END IF
    367                 ELSE ! ice present
    368                     tsurf_new(i)=t_freeze
    369                     IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
    370                         ! quantity of new ice formed over open ocean
    371                         e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
    372                                  /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
    373                         ! new ice height and fraction
    374                         h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
    375                         dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
    376                         h_new=MIN(e_freeze/dfsic,h_ice_max)
    377                         ! update tslab to freezing over open ocean only
    378                         tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki))
    379                         ! update sea ice
    380                         seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) &
    381                                    /(dfsic+fsic(ki))
    382                         fsic(ki)=fsic(ki)+dfsic
    383                         ! update snow?
    384                     END IF !freezing
    385                 END IF ! sea ice present
    386             END DO
    387         END SELECT
    388         bils_cum(:)=0.0! clear cumulated fluxes
     510          END IF
     511          IF (slab_ekman.GT.0) THEN
     512            ! dt_ekman_glo saved in W/m2
     513            dt_ekman_glo(:,:)=dt_ekman_glo(:,:)*sea_cap/cpl_pas
     514          END IF
     515        END IF ! is_mpi_root
     516!$OMP END MASTER
     517!$OMP BARRIER
     518        ! Send new fields back to all processes
     519        CALL Scatter(tslab_glo,tslab)
     520        IF (slab_hdiff) THEN
     521          CALL Scatter(dt_hdiff_glo,dt_hdiff)
     522        END IF
     523        IF (slab_ekman.GT.0) THEN
     524          CALL Scatter(dt_ekman_glo,dt_ekman)
     525          ! clear wind stress
     526          taux_cum(:)=0.
     527          tauy_cum(:)=0.
     528        END IF
     529      ENDIF ! transport
     530
     531! ***********************************
     532! Other heat fluxes
     533! ***********************************
     534      ! Add read QFlux
     535      tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas
     536      ! Add cumulated surface fluxes
     537      tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
     538      ! Convective adjustment if 2 layers
     539      IF ((nslay.GT.1).AND.(slab_cadj.GT.0)) THEN
     540        DO i=1,klon
     541          IF (tslab(i,2).GT.tslab(i,1)) THEN
     542            ! mean (mass-weighted) temperature
     543            t_cadj=SUM(tslab(i,:)*slabh(:))/SUM(slabh(:))
     544            tslab(i,1)=t_cadj
     545            tslab(i,2)=t_cadj
     546          END IF
     547        END DO
     548      END IF
     549! ***********************************
     550! Update surface temperature and ice
     551! ***********************************
     552      SELECT CASE(version_ocean)
     553      CASE('sicNO') ! no sea ice even below freezing !
     554          DO i=1,knon
     555              ki=knindex(i)
     556              tsurf_new(i)=tslab(ki,1)
     557          END DO
     558      CASE('sicOBS') ! "realistic" case, for prescribed sea ice
     559        ! tslab cannot be below freezing, or above it if there is sea ice
     560          DO i=1,knon
     561              ki=knindex(i)
     562              IF ((tslab(ki,1).LT.t_freeze).OR.(fsic(ki).GT.epsfra)) THEN
     563                  tslab(ki,1)=t_freeze
     564              END IF
     565              tsurf_new(i)=tslab(ki,1)
     566          END DO
     567      CASE('sicINT') ! interactive sea ice
     568          DO i=1,knon
     569              ki=knindex(i)
     570              IF (fsic(ki).LT.epsfra) THEN ! Free of ice
     571                  IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
     572                      ! quantity of new ice formed
     573                      e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat
     574                      ! new ice
     575                      tice(ki)=t_freeze
     576                      fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
     577                      IF (fsic(ki).GT.ice_frac_min) THEN
     578                          seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
     579                          tslab(ki,1)=t_freeze
     580                      ELSE
     581                          fsic(ki)=0.
     582                      END IF
     583                      tsurf_new(i)=t_freeze
     584                  ELSE
     585                      tsurf_new(i)=tslab(ki,1)
     586                  END IF
     587              ELSE ! ice present
     588                  tsurf_new(i)=t_freeze
     589                  IF (tslab(ki,1).LT.t_freeze) THEN ! create new ice
     590                      ! quantity of new ice formed over open ocean
     591                      e_freeze=(t_freeze-tslab(ki,1))/cyang*(1.-fsic(ki)) &
     592                               /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
     593                      ! new ice height and fraction
     594                      h_new=MIN(h_ice_new,seaice(ki)) ! max new height ice_new
     595                      dfsic=MIN(ice_frac_max-fsic(ki),e_freeze/h_new)
     596                      h_new=MIN(e_freeze/dfsic,h_ice_max)
     597                      ! update tslab to freezing over open ocean only
     598                      tslab(ki,1)=tslab(ki,1)*fsic(ki)+t_freeze*(1.-fsic(ki))
     599                      ! update sea ice
     600                      seaice(ki)=(h_new*dfsic+seaice(ki)*fsic(ki)) &
     601                                 /(dfsic+fsic(ki))
     602                      fsic(ki)=fsic(ki)+dfsic
     603                      ! update snow?
     604                  END IF ! tslab below freezing
     605              END IF ! sea ice present
     606          END DO
     607      END SELECT
     608      bils_cum(:)=0.0! clear cumulated fluxes
    389609    END IF ! coupling time
    390610  END SUBROUTINE ocean_slab_noice
    … …  
    530750           seaice(ki) = MAX(0.0,seaice(ki)-(evap(i)-snow_evap)*dtime)
    531751        ENDIF
    532 ! Melt / Freeze from above if Tsurf>0
     752! Melt / Freeze snow from above if Tsurf>0
    533753        IF (tsurf_new(i).GT.t_melt) THEN
    534             ! energy available for melting snow (in kg/m2 of snow)
     754            ! energy available for melting snow (in kg of melted snow /m2)
    535755            e_melt = MIN(MAX(snow(i)*(tsurf_new(i)-t_melt)*ice_cap/2. &
    536756               /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(i))
    … …  
    672892    END IF ! coupling time
    673893   
    674     !tests fraction glace
     894    !tests ice fraction
    675895    WHERE (fsic.LT.ice_frac_min)
    676896        tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang
    … …  
    691911    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
    692912    IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
     913    IF (ALLOCATED(tice)) DEALLOCATE(tice)
     914    IF (ALLOCATED(seaice)) DEALLOCATE(seaice)
    693915    IF (ALLOCATED(slab_bils)) DEALLOCATE(slab_bils)
    694916    IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
    695917    IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
    696918    IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum)
    697     IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
     919    IF (ALLOCATED(taux_cum)) DEALLOCATE(taux_cum)
     920    IF (ALLOCATED(tauy_cum)) DEALLOCATE(tauy_cum)
     921    IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman)
     922    IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff)
    698923
    699924  END SUBROUTINE ocean_slab_final

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

    r2471 r2669  
    30953095                endif
    30963096                mfois(nsrf) = mfois(nsrf) + 1
     3097                ! F. Codron sensible default values for ocean and sea ice
     3098                IF (nsrf.EQ.is_oce) THEN
     3099                    tsurf(i,nsrf) = 271.35 ! Freezing sea water
     3100                    DO k=1,nsw
     3101                      alb_dir(i,k,nsrf) = 0.06 ! typical Ocean albedo
     3102                      alb_dif(i,k,nsrf) = 0.06
     3103                    END DO
     3104                ELSE IF (nsrf.EQ.is_sic) THEN
     3105                    tsurf(i,nsrf) = 273.15 ! Melting ice
     3106                    DO k=1,nsw
     3107                      alb_dir(i,k,nsrf) = 0.3 ! thin ice
     3108                      alb_dif(i,k,nsrf) = 0.3
     3109                    END DO
     3110                END IF
     3111                ! F. Codron
    30973112             ELSE
    30983113                ! The continents have changed. The new fraction receives the mean sum of the existent fractions

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

    r2641 r2669  
    33SUBROUTINE phyetat0 (fichnom, clesphy0, tabcntr0)
    44
    5   USE dimphy, only: klon, zmasq, klev, nslay
     5  USE dimphy, only: klon, zmasq, klev
    66  USE iophy, ONLY : init_iophy_new
    77  USE ocean_cpl_mod,    ONLY : ocean_cpl_init
    … …  
    2525  USE carbon_cycle_mod, ONLY : carbon_cycle_tr, carbon_cycle_cpl, co2_send
    2626  USE indice_sol_mod, only: nbsrf, is_ter, epsfra, is_lic, is_oce, is_sic
    27   USE ocean_slab_mod, ONLY: tslab, seaice, tice, ocean_slab_init
     27  USE ocean_slab_mod, ONLY: nslay, tslab, seaice, tice, ocean_slab_init
    2828  USE time_phylmdz_mod, ONLY: init_iteration, pdtphys, itau_phy
    2929
    … …  
    443443  IF ( type_ocean == 'slab' ) THEN
    444444      CALL ocean_slab_init(dtime, pctsrf)
    445       found=phyetat0_get(nslay,tslab,"tslab","tslab",0.)
     445      IF (nslay.EQ.1) THEN
     446        found=phyetat0_get(1,tslab,"tslab01","tslab",0.)
     447        IF (.NOT. found) THEN
     448            found=phyetat0_get(1,tslab,"tslab","tslab",0.)
     449        END IF
     450      ELSE
     451          DO i=1,nslay
     452            WRITE(str2,'(i2.2)') i
     453            found=phyetat0_get(1,tslab(:,i),"tslab"//str2,"tslab",0.) 
     454          END DO
     455      END IF
    446456      IF (.NOT. found) THEN
    447457          PRINT*, "phyetat0: Le champ <tslab> est absent"
    … …  
    453463
    454464      ! Sea ice variables
    455       found=phyetat0_get(1,tice,"slab_tice","slab_tice",0.)
    456465      IF (version_ocean == 'sicINT') THEN
     466          found=phyetat0_get(1,tice,"slab_tice","slab_tice",0.)
    457467          IF (.NOT. found) THEN
    458468              PRINT*, "phyetat0: Le champ <tice> est absent"
    … …  
    460470                  tice(:)=ftsol(:,is_sic)
    461471          END IF
     472          found=phyetat0_get(1,seaice,"seaice","seaice",0.)
    462473          IF (.NOT. found) THEN
    463474              PRINT*, "phyetat0: Le champ <seaice> est absent"

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

    r2641 r2669  
    3030  USE indice_sol_mod, ONLY: nbsrf, is_oce, is_sic, is_ter, is_lic, epsfra
    3131  USE surface_data, ONLY: type_ocean, version_ocean
    32   USE ocean_slab_mod, ONLY : tslab, seaice, tice, fsic
     32  USE ocean_slab_mod, ONLY : nslay, tslab, seaice, tice, fsic
    3333  USE time_phylmdz_mod, ONLY: annee_ref, day_end, itau_phy, pdtphys
    3434
    … …  
    5858
    5959  INTEGER isoil, nsrf,isw
    60   CHARACTER (len=7) :: str7
     60  CHARACTER (len=2) :: str2
    6161  CHARACTER (len=256) :: nam, lnam
    6262  INTEGER           :: it, iiq
    … …  
    304304  ! Restart variables for Slab ocean
    305305  IF (type_ocean == 'slab') THEN
    306       CALL put_field("tslab", "Slab ocean temperature", tslab)
     306      IF (nslay.EQ.1) THEN
     307        CALL put_field("tslab", "Slab ocean temperature", tslab)
     308      ELSE
     309        DO it=1,nslay
     310          WRITE(str2,'(i2.2)') it
     311          CALL put_field("tslab"//str2, "Slab ocean temperature", tslab(:,it))
     312        END DO
     313      END IF
    307314      IF (version_ocean == 'sicINT') THEN
    308315          CALL put_field("seaice", "Slab seaice (kg/m2)", seaice)

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

    r2641 r2669  
    677677  TYPE(ctrl_out), SAVE :: o_slab_sic = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
    678678    'seaice', 'Epaisseur banquise slab', 'kg/m2', (/ ('', i=1, 9) /))
     679  TYPE(ctrl_out), SAVE :: o_slab_hdiff = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
     680    'slab_hdiff', 'Horizontal diffusion', 'W/m2', (/ ('', i=1, 9) /))
     681  TYPE(ctrl_out), SAVE :: o_slab_ekman = ctrl_out((/ 1, 1, 10, 10, 10, 10, 11, 11, 11 /), &
     682    'slab_ekman', 'Ekman heat transport', 'W/m2', (/ ('', i=1, 9) /))
    679683  TYPE(ctrl_out), SAVE :: o_ale_bl = ctrl_out((/ 1, 1, 1, 10, 10, 10, 11, 11, 11 /), &
    680684    'ale_bl', 'ALE BL', 'm2/s2', (/ ('', i=1, 9) /))

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

    r2641 r2669  
    2424    ! defined and initialised in phys_output_mod.F90
    2525
    26     USE dimphy, only: klon, klev, klevp1, nslay
     26    USE dimphy, only: klon, klev, klevp1
    2727    USE mod_phys_lmdz_para, ONLY: is_north_pole_phy,is_south_pole_phy
    2828    USE mod_grid_phy_lmdz, ONLY : nbp_lon, nbp_lat
    … …  
    8989         o_slab_qflux, o_tslab, o_slab_bils, &
    9090         o_slab_bilg, o_slab_sic, o_slab_tice, &
     91         o_slab_hdiff, o_slab_ekman, &
    9192         o_weakinv, o_dthmin, o_cldtau, &
    9293         o_cldemi, o_pr_con_l, o_pr_con_i, &
    … …  
    291292 
    292293
    293     USE ocean_slab_mod, only: tslab, slab_bils, slab_bilg, tice, seaice
     294    USE ocean_slab_mod, only: nslay, tslab, slab_bils, slab_bilg, tice, &
     295        seaice, slab_ekman,slab_hdiff, dt_ekman, dt_hdiff
    294296    USE pbl_surface_mod, only: snow
    295297    USE indice_sol_mod, only: nbsrf
    … …  
    312314    USE phys_cal_mod, only : mth_len
    313315
     316#ifdef CPP_RRTM
     317    USE YOESW, ONLY : RSUN
     318#endif
    314319
    315320    IMPLICIT NONE
    … …  
    355360#endif
    356361    REAL, PARAMETER :: un_jour=86400.
     362    INTEGER ISW
     363    CHARACTER*1 ch1
    357364
    358365    ! On calcul le nouveau tau:
    … …  
    403410       CALL histwrite_phy(o_contfracOR, pctsrf(:,is_ter))
    404411       CALL histwrite_phy(o_aireTER, paire_ter)
     412!
     413#ifdef CPP_XIOS
     414       CALL histwrite_phy("R_ecc",R_ecc)
     415       CALL histwrite_phy("R_peri",R_peri)
     416       CALL histwrite_phy("R_incl",R_incl)
     417       CALL histwrite_phy("solaire",solaire)
     418!
     419#ifdef CPP_RRTM
     420      IF (iflag_rrtm.EQ.1) THEN
     421        DO ISW=1, NSW
     422          WRITE(ch1,'(i1)') ISW
     423!         zx_tmp_0d=RSUN(ISW)
     424!         CALL histwrite_phy("rsun"//ch1,zx_tmp_0d)
     425          CALL histwrite_phy("rsun"//ch1,RSUN(ISW))
     426        ENDDO
     427      ENDIF
     428#endif
     429!
     430       CALL histwrite_phy("co2_ppm",co2_ppm)
     431       CALL histwrite_phy("CH4_ppb",CH4_ppb)
     432       CALL histwrite_phy("N2O_ppb",N2O_ppb)
     433       CALL histwrite_phy("CFC11_ppt",CFC11_ppt)
     434       CALL histwrite_phy("CFC12_ppt",CFC12_ppt)
     435!
     436#endif
     437
    405438!!! Champs 2D !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    406439! Simulateur AIRS
    … …  
    951984              CALL histwrite_phy(o_tslab, zx_tmp_fi2d)
    952985          ELSE
    953               CALL histwrite_phy(o_tslab, tslab)
     986              CALL histwrite_phy(o_tslab, tslab(:,1:nslay))
    954987          END IF
    955988          IF (version_ocean=='sicINT') THEN
    … …  
    957990              CALL histwrite_phy(o_slab_tice, tice)
    958991              CALL histwrite_phy(o_slab_sic, seaice)
     992          END IF
     993          IF (slab_hdiff) THEN
     994            IF (nslay.EQ.1) THEN
     995                zx_tmp_fi2d(:)=dt_hdiff(:,1)
     996                CALL histwrite_phy(o_slab_hdiff, zx_tmp_fi2d)
     997            ELSE
     998                CALL histwrite_phy(o_slab_hdiff, dt_hdiff(:,1:nslay))
     999            END IF
     1000          END IF
     1001          IF (slab_ekman.GT.0) THEN
     1002            IF (nslay.EQ.1) THEN
     1003                zx_tmp_fi2d(:)=dt_ekman(:,1)
     1004                CALL histwrite_phy(o_slab_ekman, zx_tmp_fi2d)
     1005            ELSE
     1006                CALL histwrite_phy(o_slab_ekman, dt_ekman(:,1:nslay))
     1007            END IF
    9591008          END IF
    9601009       ENDIF !type_ocean == force/slab

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

    r2641 r2669  
    233233    USE ioipsl_getin_p_mod, ONLY : getin_p
    234234
     235#ifndef CPP_XIOS
    235236    USE paramLMDZ_phy_mod
     237#endif
    236238
    237239    USE cmp_seri_mod
    … …  
    582584    LOGICAL,SAVE :: ok_adjwk=.FALSE.
    583585    !$OMP THREADPRIVATE(ok_adjwk)
    584     REAL,SAVE :: oliqmax=999.
    585     !$OMP THREADPRIVATE(oliqmax)
     586    REAL,SAVE :: oliqmax=999.,oicemax=999.
     587    !$OMP THREADPRIVATE(oliqmax,oicemax)
    586588    REAL, SAVE :: alp_offset
    587589    !$OMP THREADPRIVATE(alp_offset)
    … …  
    983985    INTEGER, SAVE :: flag_aerosol
    984986    !$OMP THREADPRIVATE(flag_aerosol)
     987    LOGICAL, SAVE :: flag_bc_internal_mixture
     988    !$OMP THREADPRIVATE(flag_bc_internal_mixture)
    985989    LOGICAL, SAVE :: new_aod
    986990    !$OMP THREADPRIVATE(new_aod)
    … …  
    10111015    ! edges of pressure intervals for ozone climatologies, in Pa, in strictly
    10121016    ! ascending order
    1013 
    1014     integer, save:: co3i = 0
    1015     !     time index in NetCDF file of current ozone fields
    1016     !$OMP THREADPRIVATE(co3i)
    10171017
    10181018    integer ro3i
    … …  
    11401140            ok_ade, ok_aie, ok_cdnc, aerosol_couple,  &
    11411141            flag_aerosol, flag_aerosol_strat, new_aod, &
    1142             bl95_b0, bl95_b1, &
     1142            flag_bc_internal_mixture, bl95_b0, bl95_b1, &
    11431143                                ! nv flags pour la convection et les
    11441144                                ! poches froides
    … …  
    11971197       CALL getin_p('ok_adjwk',ok_adjwk)
    11981198       CALL getin_p('oliqmax',oliqmax)
     1199       CALL getin_p('oicemax',oicemax)
    11991200       CALL getin_p('ratqsp0',ratqsp0)
    12001201       CALL getin_p('ratqsdp',ratqsdp)
    … …  
    15511552       WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
    15521553
     1554#ifndef CPP_XIOS
    15531555       CALL ini_paramLMDZ_phy(dtime,nid_ctesGCM)
     1556#endif
    15541557
    15551558#endif
    … …  
    19061909    ! Prescrire l'ozone et calculer l'albedo sur l'ocean.
    19071910    !
    1908     if (read_climoz >= 1) then
    1909        ! Ozone from a file
    1910        ! Update required ozone index:
    1911        ro3i = int((days_elapsed + jh_cur - jh_1jan) / year_len * 360.) + 1
    1912        if (ro3i == 361) ro3i = 360
    1913        ! (This should never occur, except perhaps because of roundup
    1914        ! error. See documentation.)
    1915        if (ro3i /= co3i) then
    1916           ! Update ozone field:
    1917           if (read_climoz == 1) then
    1918              call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
    1919                   press_in_edg=press_climoz, paprs=paprs, v3=wo)
    1920           else
    1921              ! read_climoz == 2
    1922              call regr_pr_av(ncid_climoz, (/"tro3         ", &
    1923                   "tro3_daylight"/), julien=ro3i, press_in_edg=press_climoz, &
    1924                   paprs=paprs, v3=wo)
    1925           end if
    1926           ! Convert from mole fraction of ozone to column density of ozone in a
    1927           ! cell, in kDU:
    1928           forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
    1929                * zmasse / dobson_u / 1e3
    1930           ! (By regridding ozone values for LMDZ only once every 360th of
    1931           ! year, we have already neglected the variation of pressure in one
    1932           ! 360th of year. So do not recompute "wo" at each time step even if
    1933           ! "zmasse" changes a little.)
    1934           co3i = ro3i
    1935        end if
    1936     ELSEIF (MOD(itap-1,lmt_pas) == 0) THEN
     1911    ! Update ozone if day change
     1912    IF (MOD(itap-1,lmt_pas) == 0) THEN
     1913      IF (read_climoz == 0) THEN
    19371914       ! Once per day, update ozone from Royer:
    1938 
    1939        IF (solarlong0<-999.) then
    1940           ! Generic case with evolvoing season
    1941           zzz=real(days_elapsed+1)
    1942        ELSE IF (abs(solarlong0-1000.)<1.e-4) then
    1943           ! Particular case with annual mean insolation
    1944           zzz=real(90) ! could be revisited
    1945           IF (read_climoz/=-1) THEN
    1946              abort_message ='read_climoz=-1 is recommended when ' &
    1947                   // 'solarlong0=1000.'
    1948              CALL abort_physic (modname,abort_message,1)
    1949           ENDIF
    1950        ELSE
     1915        IF (solarlong0<-999.) then
     1916           ! Generic case with evolvoing season
     1917           zzz=real(days_elapsed+1)
     1918        ELSE IF (abs(solarlong0-1000.)<1.e-4) then
     1919           ! Particular case with annual mean insolation
     1920           zzz=real(90) ! could be revisited
     1921           IF (read_climoz/=-1) THEN
     1922              abort_message ='read_climoz=-1 is recommended when ' &
     1923                   // 'solarlong0=1000.'
     1924              CALL abort_physic (modname,abort_message,1)
     1925           ENDIF
     1926        ELSE
    19511927          ! Case where the season is imposed with solarlong0
    19521928          zzz=real(90) ! could be revisited
    1953        ENDIF
    1954        wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
     1929        ENDIF
     1930
     1931        wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
     1932      ELSE
     1933        ro3i = int((days_elapsed + jh_cur - jh_1jan) / year_len * 360.) + 1   
     1934        if (ro3i == 361) ro3i = 360
     1935        if (read_climoz == 1) then
     1936           call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
     1937                press_in_edg=press_climoz, paprs=paprs, v3=wo)
     1938        else
     1939           ! read_climoz == 2
     1940           call regr_pr_av(ncid_climoz, (/"tro3         ", &
     1941                "tro3_daylight"/), julien=ro3i, press_in_edg=press_climoz, &
     1942                paprs=paprs, v3=wo)
     1943        end if
     1944        ! Convert from mole fraction of ozone to column density of ozone in a
     1945        ! cell, in kDU:
     1946        forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
     1947             * zmasse / dobson_u / 1e3
     1948        ! (By regridding ozone values for LMDZ only once every 360th of
     1949        ! year, we have already neglected the variation of pressure in one
     1950        ! 360th of year. So do not recompute "wo" at each time step even if
     1951        ! "zmasse" changes a little.)
     1952     
     1953      ENDIF
    19551954    ENDIF
    19561955    !
    … …  
    30343033         'lsc',abortphy,flag_inhib_tend)
    30353034    rain_num(:)=0.
    3036     DO k = 1, klev  ! re-evaporation de l'eau liquide nuageuse
     3035    DO k = 1, klev
    30373036       DO i = 1, klon
    30383037          IF (ql_seri(i,k)>oliqmax) THEN
    … …  
    30423041       ENDDO
    30433042    ENDDO
     3043    IF (nqo==3) THEN
     3044    DO k = 1, klev
     3045       DO i = 1, klon
     3046          IF (qs_seri(i,k)>oicemax) THEN
     3047             rain_num(i)=rain_num(i)+(qs_seri(i,k)-oicemax)*zmasse(i,k)/pdtphys
     3048             qs_seri(i,k)=oicemax
     3049          ENDIF
     3050       ENDDO
     3051    ENDDO
     3052    ENDIF
    30443053
    30453054    !---------------------------------------------------------------------------
    … …  
    33883397                   !
    33893398                   CALL readaerosol_optic_rrtm( debut, aerosol_couple, &
    3390                         new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
     3399                        new_aod, flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
    33913400                        pdtphys, pplay, paprs, t_seri, rhcl, presnivs,  &
    33923401                        tr_seri, mass_solu_aero, mass_solu_aero_pi,  &
    … …  
    45174526
    45184527
     4528#ifndef CPP_XIOS
    45194529    CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
     4530#endif
    45204531
    45214532#endif

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

    r2408 r2669  
    166166  m_allaer_pi(:,:,id_CSSO4M_phy)  = 0.                   ! CSSO4M pre-ind
    167167  m_allaer_pi(:,:,id_SSSSM_phy)   = sssupco_pi(:,:)      ! SSSSM pre-ind
    168   m_allaer_pi(:,:,id_ASSSM_phy)   = sscoarse_pi(:,:)     ! CSSSM pre-ind
    169   m_allaer_pi(:,:,id_CIDUSTM_phy) = ssacu_pi(:,:)        ! ASSSM pre-ind
    170   m_allaer_pi(:,:,id_AIBCM_phy)  = cidust_pi(:,:)       ! CIDUSTM pre-ind
     168  m_allaer_pi(:,:,id_CSSSM_phy)   = sscoarse_pi(:,:)     ! CSSSM pre-ind
     169  m_allaer_pi(:,:,id_ASSSM_phy)  = ssacu_pi(:,:)        ! ASSSM pre-ind
     170  m_allaer_pi(:,:,id_CIDUSTM_phy) = cidust_pi(:,:)       ! CIDUSTM pre-ind
    171171  m_allaer_pi(:,:,id_AIBCM_phy)   = bcins_pi(:,:)        ! AIBCM pre-ind
    172172  m_allaer_pi(:,:,id_AIPOMM_phy)  = pomins_pi(:,:)       ! AIPOMM pre-ind

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

    r2641 r2669  
    33!
    44
    5 SUBROUTINE AEROPT_5WV_RRTM(&
    6    pdel, m_allaer, delt, &
    7    RHcl, ai, flag_aerosol, &
    8    pplay, t_seri, &
     5SUBROUTINE AEROPT_5WV_RRTM(  &
     6   pdel, m_allaer,           &
     7   RHcl, ai, flag_aerosol,   &
     8   flag_bc_internal_mixture, &
     9   pplay, t_seri,            &
    910   tausum, tau )
    1011
    … …  
    5455  !
    5556  REAL, DIMENSION(klon,klev), INTENT(in)   :: pdel
    56   REAL, INTENT(in)                         :: delt
    5757  REAL, DIMENSION(klon,klev,naero_tot), INTENT(in) :: m_allaer
    5858  REAL, DIMENSION(klon,klev), INTENT(in)   :: RHcl     ! humidite relative ciel clair
    5959  INTEGER,INTENT(in)                       :: flag_aerosol
     60  LOGICAL,INTENT(in)                       :: flag_bc_internal_mixture
    6061  REAL, DIMENSION(klon,klev), INTENT(in)   :: pplay
    6162  REAL, DIMENSION(klon,klev), INTENT(in)   :: t_seri
    … …  
    106107
    107108  !
     109  ! BC internal mixture
     110  !
     111  INTEGER, PARAMETER ::  nbclassbc = 5  ! Added by Rong Wang/OB for the 5 fractions
     112                                       ! of BC in the soluble mode:
     113                                       ! bc_content/0.001, 0.01, 0.02, 0.05, ! 0.1/
     114  ! for Maxwell-Garnet internal mixture
     115  ! Detailed theory can be found in R. Wang Estimation of global black carbon ! direct
     116  ! radiative forcing and its uncertainty constrained by observations. J.
     117  ! Geophys. Res. Atmos. Added by R. Wang and OB
     118  REAL :: alpha_MG_5wv(nbre_RH,las,nbclassbc)
     119
     120  !
    108121  ! Proprietes optiques
    109122  !
    110   REAL :: fact_RH(nbre_RH)
    111   INTEGER :: n
     123  REAL :: fact_RH(nbre_RH), BC_massfra
     124  INTEGER :: n, classbc
    112125
    113126! From here on we look at the optical parameters at 5 wavelengths: 
    … …  
    117130 DATA alpha_aers_5wv/ &
    118131                                ! bc soluble
    119        7.930,7.930,7.930,7.930,7.930,7.930,     &
    120        7.930,7.930,10.893,12.618,14.550,16.613, &
    121        7.658,7.658,7.658,7.658,7.658,7.658,     &
    122        7.658,7.658,10.351,11.879,13.642,15.510, &
    123        7.195,7.195,7.195,7.195,7.195,7.195,     &
    124        7.195,7.195,9.551,10.847,12.381,13.994,  &
    125        6.736,6.736,6.736,6.736,6.736,6.736,     &
    126        6.736,6.736,8.818,9.938,11.283,12.687,   &
    127        6.277,6.277,6.277,6.277,6.277,6.277,     &
    128        6.277,6.277,8.123,9.094,10.275,11.501,   &
     132       7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, &
     133       7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, &
     134       7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994,  &
     135       6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687,   &
     136       6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501,   &
    129137                                ! pom soluble
    130        6.676,6.676,6.676,6.676,6.710,6.934,   &
    131        7.141,7.569,8.034,8.529,9.456,10.511,  &
    132        5.109,5.109,5.109,5.109,5.189,5.535,   &
    133        5.960,6.852,8.008,9.712,12.897,19.676, &
    134        3.718,3.718,3.718,3.718,3.779,4.042,   &
    135        4.364,5.052,5.956,7.314,9.896,15.688,  &
    136        2.849,2.849,2.849,2.849,2.897,3.107,   &
    137        3.365,3.916,4.649,5.760,7.900,12.863,  &
    138        2.229,2.229,2.229,2.229,2.268,2.437,   &
    139        2.645,3.095,3.692,4.608,6.391,10.633,  &
     138       6.676,6.676,6.676,6.676,6.710,6.934,7.141,7.569,8.034,8.529,9.456,10.511,  &
     139       5.109,5.109,5.109,5.109,5.189,5.535,5.960,6.852,8.008,9.712,12.897,19.676, &
     140       3.718,3.718,3.718,3.718,3.779,4.042,4.364,5.052,5.956,7.314,9.896,15.688,  &
     141       2.849,2.849,2.849,2.849,2.897,3.107,3.365,3.916,4.649,5.760,7.900,12.863,  &
     142       2.229,2.229,2.229,2.229,2.268,2.437,2.645,3.095,3.692,4.608,6.391,10.633,  &
    140143                                ! Sulfate (Accumulation)
    141        5.751,6.215,6.690,7.024,7.599,8.195,      &
    142        9.156,10.355,12.660,14.823,18.908,24.508, &
    143        4.320,4.675,5.052,5.375,5.787,6.274,      &
    144        7.066,8.083,10.088,12.003,15.697,21.133,  &
    145        3.079,3.351,3.639,3.886,4.205,4.584,      &
    146        5.206,6.019,7.648,9.234,12.391,17.220,    &
    147        2.336,2.552,2.781,2.979,3.236,3.540,      &
    148        4.046,4.711,6.056,7.388,10.093,14.313,    &
    149        1.777,1.949,2.134,2.292,2.503,2.751,      &
    150        3.166,3.712,4.828,5.949,8.264,11.922,     &
     144       5.751,6.215,6.690,7.024,7.599,8.195,9.156,10.355,12.660,14.823,18.908,24.508, &
     145       4.320,4.675,5.052,5.375,5.787,6.274,7.066,8.083,10.088,12.003,15.697,21.133,  &
     146       3.079,3.351,3.639,3.886,4.205,4.584,5.206,6.019,7.648,9.234,12.391,17.220,    &
     147       2.336,2.552,2.781,2.979,3.236,3.540,4.046,4.711,6.056,7.388,10.093,14.313,    &
     148       1.777,1.949,2.134,2.292,2.503,2.751,3.166,3.712,4.828,5.949,8.264,11.922,     &
    151149                                ! Sulfate (Coarse)
    152        5.751,6.215,6.690,7.024,7.599,8.195,      &
    153        9.156,10.355,12.660,14.823,18.908,24.508, &
    154        4.320,4.675,5.052,5.375,5.787,6.274,      &
    155        7.066,8.083,10.088,12.003,15.697,21.133,  &
    156        3.079,3.351,3.639,3.886,4.205,4.584,      &
    157        5.206,6.019,7.648,9.234,12.391,17.220,    &
    158        2.336,2.552,2.781,2.979,3.236,3.540,      &
    159        4.046,4.711,6.056,7.388,10.093,14.313,    &
    160        1.777,1.949,2.134,2.292,2.503,2.751,      &
    161        3.166,3.712,4.828,5.949,8.264,11.922,     &
     150       5.751,6.215,6.690,7.024,7.599,8.195,9.156,10.355,12.660,14.823,18.908,24.508, &
     151       4.320,4.675,5.052,5.375,5.787,6.274,7.066,8.083,10.088,12.003,15.697,21.133,  &
     152       3.079,3.351,3.639,3.886,4.205,4.584,5.206,6.019,7.648,9.234,12.391,17.220,    &
     153       2.336,2.552,2.781,2.979,3.236,3.540,4.046,4.711,6.056,7.388,10.093,14.313,    &
     154       1.777,1.949,2.134,2.292,2.503,2.751,3.166,3.712,4.828,5.949,8.264,11.922,     &
    162155                           ! seasalt seasalt Super Coarse Soluble (SS)
    163         0.218, 0.272, 0.293, 0.316, 0.343, 0.380, &
    164         0.429, 0.501, 0.636, 0.755, 0.967, 1.495, &
    165         0.221, 0.275, 0.297, 0.320, 0.348, 0.383, &
    166         0.432, 0.509, 0.640, 0.759, 0.972, 1.510, &
    167         0.224, 0.279, 0.301, 0.324, 0.352, 0.388, &
    168         0.438, 0.514, 0.647, 0.768, 0.985, 1.514, &
    169         0.227, 0.282, 0.303, 0.327, 0.356, 0.392, &
    170         0.441, 0.518, 0.652, 0.770, 0.987, 1.529, &
    171         0.230, 0.285, 0.306, 0.330, 0.359, 0.396, &
    172         0.446, 0.522, 0.656, 0.777, 0.993, 1.539, &
     156        0.218, 0.272, 0.293, 0.316, 0.343, 0.380, 0.429, 0.501, 0.636, 0.755, 0.967, 1.495, &
     157        0.221, 0.275, 0.297, 0.320, 0.348, 0.383, 0.432, 0.509, 0.640, 0.759, 0.972, 1.510, &
     158        0.224, 0.279, 0.301, 0.324, 0.352, 0.388, 0.438, 0.514, 0.647, 0.768, 0.985, 1.514, &
     159        0.227, 0.282, 0.303, 0.327, 0.356, 0.392, 0.441, 0.518, 0.652, 0.770, 0.987, 1.529, &
     160        0.230, 0.285, 0.306, 0.330, 0.359, 0.396, 0.446, 0.522, 0.656, 0.777, 0.993, 1.539, &
    173161                           ! seasalt seasalt Coarse Soluble (CS)     
    174         0.578, 0.706, 0.756, 0.809, 0.876, 0.964, &
    175         1.081, 1.256, 1.577, 1.858, 2.366, 3.613, &
    176         0.598, 0.725, 0.779, 0.833, 0.898, 0.990, &
    177         1.109, 1.290, 1.609, 1.889, 2.398, 3.682, &
    178         0.619, 0.750, 0.802, 0.857, 0.927, 1.022, &
    179         1.141, 1.328, 1.648, 1.939, 2.455, 3.729, &
    180         0.633, 0.767, 0.820, 0.879, 0.948, 1.044, &
    181         1.167, 1.353, 1.683, 1.969, 2.491, 3.785, &
    182         0.648, 0.785, 0.838, 0.896, 0.967, 1.066, &
    183         1.192, 1.381, 1.714, 2.006, 2.531, 3.836, &
     162        0.578, 0.706, 0.756, 0.809, 0.876, 0.964, 1.081, 1.256, 1.577, 1.858, 2.366, 3.613, &
     163        0.598, 0.725, 0.779, 0.833, 0.898, 0.990, 1.109, 1.290, 1.609, 1.889, 2.398, 3.682, &
     164        0.619, 0.750, 0.802, 0.857, 0.927, 1.022, 1.141, 1.328, 1.648, 1.939, 2.455, 3.729, &
     165        0.633, 0.767, 0.820, 0.879, 0.948, 1.044, 1.167, 1.353, 1.683, 1.969, 2.491, 3.785, &
     166        0.648, 0.785, 0.838, 0.896, 0.967, 1.066, 1.192, 1.381, 1.714, 2.006, 2.531, 3.836, &
    184167                           ! seasalt seasalt Accumulation Soluble (AS)
    185         4.432, 5.899, 6.505, 7.166, 7.964, 7.962, &
    186         9.232,11.257,14.979,18.337,24.223,37.811, &
    187         3.298, 4.569, 5.110, 5.709, 6.446, 6.268, &
    188         7.396, 9.246,12.787,16.113,22.197,37.136, &
    189         2.340, 3.358, 3.803, 4.303, 4.928, 4.696, &
    190         5.629, 7.198,10.308,13.342,19.120,34.296, &
    191         1.789, 2.626, 2.999, 3.422, 3.955, 3.730, &
    192         4.519, 5.864, 8.593,11.319,16.653,31.331, &
    193         1.359, 2.037, 2.343, 2.693, 3.139, 2.940, &
    194         3.596, 4.729, 7.076, 9.469,14.266,28.043 /
     168        4.432, 5.899, 6.505, 7.166, 7.964, 7.962, 9.232,11.257,14.979,18.337,24.223,37.811, &
     169        3.298, 4.569, 5.110, 5.709, 6.446, 6.268, 7.396, 9.246,12.787,16.113,22.197,37.136, &
     170        2.340, 3.358, 3.803, 4.303, 4.928, 4.696, 5.629, 7.198,10.308,13.342,19.120,34.296, &
     171        1.789, 2.626, 2.999, 3.422, 3.955, 3.730, 4.519, 5.864, 8.593,11.319,16.653,31.331, &
     172        1.359, 2.037, 2.343, 2.693, 3.139, 2.940, 3.596, 4.729, 7.076, 9.469,14.266,28.043 /
    195173
    196174  DATA alpha_aeri_5wv/ &
    … …  
    201179                                 ! pom insoluble
    202180        5.042, 3.101, 1.890, 1.294, 0.934/
     181
     182  ! internal mixture
     183  DATA alpha_MG_5wv/ &
     184     ! bc content = 0.001
     185       7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, &
     186       7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, &
     187       7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994,  &
     188       6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687,   &
     189       6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501,   &
     190     ! bc content = 0.01
     191       7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, &
     192       7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, &
     193       7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994,  &
     194       6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687,   &
     195       6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501,   &
     196     ! bc content = 0.02
     197       7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, &
     198       7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, &
     199       7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994,  &
     200       6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687,   &
     201       6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501,   &
     202     ! bc content = 0.05
     203       7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, &
     204       7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, &
     205       7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994,  &
     206       6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687,   &
     207       6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501,   &
     208     ! bc content = 0.10
     209       7.930,7.930,7.930,7.930,7.930,7.930,7.930,7.930,10.893,12.618,14.550,16.613, &
     210       7.658,7.658,7.658,7.658,7.658,7.658,7.658,7.658,10.351,11.879,13.642,15.510, &
     211       7.195,7.195,7.195,7.195,7.195,7.195,7.195,7.195,9.551,10.847,12.381,13.994,  &
     212       6.736,6.736,6.736,6.736,6.736,6.736,6.736,6.736,8.818,9.938,11.283,12.687,   &
     213       6.277,6.277,6.277,6.277,6.277,6.277,6.277,6.277,8.123,9.094,10.275,11.501 /
     214
    203215  !
    204216  ! Initialisations
    … …  
    251263     aerosol_name(8) = id_ASSSM_phy
    252264     aerosol_name(9) = id_CIDUSTM_phy
    253      aerosol_name(10) = id_CSSO4M_phy
     265     aerosol_name(10)= id_CSSO4M_phy
    254266  ENDIF
    255267
    … …  
    289301        soluble=.TRUE.
    290302        spsol=3
    291         fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD
     303        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for AOD
    292304    ELSEIF (aerosol_name(m).EQ.id_CSSO4M_phy) THEN
    293305        soluble=.TRUE.
    294306        spsol=4
    295         fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD
     307        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for AOD
    296308    ELSEIF (aerosol_name(m).EQ.id_SSSSM_phy) THEN
    297309        soluble=.TRUE.
    … …  
    325337      IF (soluble) THEN            ! For soluble aerosol
    326338
     339        !--treat special case of soluble BC internal mixture
     340        IF (spsol.EQ.1 .AND. flag_bc_internal_mixture) THEN
     341
     342          DO k=1, klev
     343            DO i=1, klon
     344
     345             BC_massfra = m_allaer(i,k,id_ASBCM_phy)/(m_allaer(i,k,id_ASBCM_phy)+m_allaer(i,k,id_ASSO4M_phy))
     346
     347             IF (BC_massfra.GE.0.10) THEN
     348               classbc = 5
     349             ELSEIF  (BC_massfra.GE.0.05) THEN
     350               classbc = 4
     351             ELSEIF  (BC_massfra.GE.0.02) THEN
     352               classbc = 3
     353             ELSEIF  (BC_massfra.GE.0.01) THEN
     354               classbc = 2
     355             ELSE
     356               classbc = 1
     357             ENDIF
     358
     359              tau_ae5wv_int = alpha_MG_5wv(RH_num(i,k),la,classbc)+DELTA(i,k)* &
     360                             (alpha_MG_5wv(RH_num(i,k)+1,la,classbc) - &
     361                              alpha_MG_5wv(RH_num(i,k),la,classbc))
     362              tau(i,k,la,aerindex) = m_allaer(i,k,aerindex)/1.e6*zdh(i,k)*tau_ae5wv_int*fac
     363              tausum(i,la,aerindex)=tausum(i,la,aerindex)+tau(i,k,la,aerindex)
     364            ENDDO
     365          ENDDO
     366
     367        !--other cases of soluble aerosols
     368        ELSE
     369
    327370          DO k=1, klev
    328371            DO i=1, klon
    … …  
    334377            ENDDO
    335378          ENDDO
     379
     380        ENDIF
    336381 
    337       ELSE                         ! For insoluble aerosol
     382      ! cases of insoluble aerosol
     383      ELSE                         
    338384
    339385        DO k=1, klev

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

    r2641 r2669  
    33!
    44SUBROUTINE AEROPT_6BANDS_RRTM ( &
    5      pdel, m_allaer, delt, RHcl, &
     5     pdel, m_allaer, RHcl, &
    66     tau_allaer, piz_allaer, &
    77     cg_allaer, m_allaer_pi, &
    8      flag_aerosol, zrho )
     8     flag_aerosol, flag_bc_internal_mixture, zrho )
    99
    1010  USE dimphy
    … …  
    2828  !
    2929  REAL, DIMENSION(klon,klev),     INTENT(in)  :: pdel
    30   REAL,                           INTENT(in)  :: delt
    3130  REAL, DIMENSION(klon,klev,naero_tot),   INTENT(in)  :: m_allaer
    3231  REAL, DIMENSION(klon,klev,naero_tot),   INTENT(in)  :: m_allaer_pi
    3332  REAL, DIMENSION(klon,klev),     INTENT(in)  :: RHcl       ! humidite relative ciel clair
    3433  INTEGER,                        INTENT(in)  :: flag_aerosol
     34  LOGICAL,                        INTENT(in)  :: flag_bc_internal_mixture
    3535  REAL, DIMENSION(klon,klev),     INTENT(in)  :: zrho
    3636  !
    … …  
    7171
    7272  REAL, DIMENSION(klon,klev,naero_tot,nbands_sw_rrtm) ::  tau_ae
    73   REAL, DIMENSION(klon,klev,naero_tot,nbands_sw_rrtm) ::  tau_ae_pi
    7473  REAL, DIMENSION(klon,klev,naero_tot,nbands_sw_rrtm) ::  piz_ae
    7574  REAL, DIMENSION(klon,klev,naero_tot,nbands_sw_rrtm) ::  cg_ae
     75
     76  REAL, DIMENSION(klon,klev,naero_tot,nbands_sw_rrtm) ::  tau_ae_pi
     77  REAL, DIMENSION(klon,klev,id_ASBCM_phy:id_ASBCM_phy,nbands_sw_rrtm) :: piz_ae_pi
     78  REAL, DIMENSION(klon,klev,id_ASBCM_phy:id_ASBCM_phy,nbands_sw_rrtm) :: cg_ae_pi
    7679  !
    7780  ! Proprietes optiques
    … …  
    8386  REAL:: piz_aers_6bands(nbre_RH,nbands_sw_rrtm,naero_soluble)     !-- unit
    8487  REAL:: piz_aeri_6bands(nbands_sw_rrtm,naero_insoluble)        !-- unit
    85 
    86   INTEGER :: id
    87   REAL :: tmp_var, tmp_var_pi
     88  !
     89  ! BC internal mixture
     90  !
     91  INTEGER, PARAMETER ::  nbclassbc = 5  ! Added by Rong Wang/OB for the 5 fractions
     92                                       ! of BC in the soluble mode:
     93                                       ! bc_content/0.001, 0.01, 0.02, 0.05, 0.1/
     94  ! for Maxwell-Garnet internal mixture
     95  ! Detailed theory can be found in R. Wang Estimation of global black carbon direct
     96  ! radiative forcing and its uncertainty constrained by observations. J.
     97  ! Geophys. Res. Atmos. Added by R. Wang and OB
     98  REAL :: alpha_MG_6bands(nbre_RH,nbands_sw_rrtm,nbclassbc)
     99  REAL :: cg_MG_6bands(nbre_RH,nbands_sw_rrtm,nbclassbc)
     100  REAL :: piz_MG_6bands(nbre_RH,nbands_sw_rrtm,nbclassbc)
     101  !
     102  INTEGER :: id, classbc, classbc_pi
     103  REAL :: tmp_var, tmp_var_pi, BC_massfra, BC_massfra_pi
     104
     105  !
     106  REAL, PARAMETER :: tau_min = 1.e-15
     107!  REAL, PARAMETER :: tau_min = 1.e-7
    88108
    89109!***************************************************************************
    … …  
    268288  0.973, 0.973, 0.972, 0.940, 0.816, 0.663 /
    269289
     290! Added by R. Wang (July 31 2016)
     291! properties for BC assuming Maxwell-Garnett rule and internal mixture
     292
     293  DATA alpha_MG_6bands/ &
     294     ! bc content = 0.001
     295  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
     296  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
     297  5.900, 5.900, 5.900, 5.900, 5.900, 6.502, 7.151, 8.496, 9.818, 9.965,10.124,11.564, &
     298  4.284, 4.284, 4.284, 4.284, 4.284, 4.721, 5.193, 6.169, 7.129, 7.236, 7.352, 8.397, &
     299  2.163, 2.163, 2.163, 2.163, 2.163, 2.384, 2.622, 3.115, 3.600, 3.654, 3.712, 4.240, &
     300  0.966, 0.966, 0.966, 0.966, 0.966, 1.065, 1.171, 1.392, 1.608, 1.632, 1.658, 1.894, &
     301     ! bc content = 0.01
     302  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
     303  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
     304  5.900, 5.900, 5.900, 5.900, 5.900, 6.502, 7.151, 8.496, 9.818, 9.965,10.124,11.564, &
     305  4.284, 4.284, 4.284, 4.284, 4.284, 4.721, 5.193, 6.169, 7.129, 7.236, 7.352, 8.397, &
     306  2.163, 2.163, 2.163, 2.163, 2.163, 2.384, 2.622, 3.115, 3.600, 3.654, 3.712, 4.240, &
     307  0.966, 0.966, 0.966, 0.966, 0.966, 1.065, 1.171, 1.392, 1.608, 1.632, 1.658, 1.894, &
     308     ! bc content = 0.02
     309  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
     310  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
     311  5.900, 5.900, 5.900, 5.900, 5.900, 6.502, 7.151, 8.496, 9.818, 9.965,10.124,11.564, &
     312  4.284, 4.284, 4.284, 4.284, 4.284, 4.721, 5.193, 6.169, 7.129, 7.236, 7.352, 8.397, &
     313  2.163, 2.163, 2.163, 2.163, 2.163, 2.384, 2.622, 3.115, 3.600, 3.654, 3.712, 4.240, &
     314  0.966, 0.966, 0.966, 0.966, 0.966, 1.065, 1.171, 1.392, 1.608, 1.632, 1.658, 1.894, &
     315     ! bc content = 0.05
     316  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
     317  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
     318  5.900, 5.900, 5.900, 5.900, 5.900, 6.502, 7.151, 8.496, 9.818, 9.965,10.124,11.564, &
     319  4.284, 4.284, 4.284, 4.284, 4.284, 4.721, 5.193, 6.169, 7.129, 7.236, 7.352, 8.397, &
     320  2.163, 2.163, 2.163, 2.163, 2.163, 2.384, 2.622, 3.115, 3.600, 3.654, 3.712, 4.240, &
     321  0.966, 0.966, 0.966, 0.966, 0.966, 1.065, 1.171, 1.392, 1.608, 1.632, 1.658, 1.894, &
     322     ! bc content = 0.10
     323  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
     324  6.497, 6.497, 6.497, 6.497, 6.497, 7.160, 7.875, 9.356,10.811,10.974,11.149,12.734, &
     325  5.900, 5.900, 5.900, 5.900, 5.900, 6.502, 7.151, 8.496, 9.818, 9.965,10.124,11.564, &
     326  4.284, 4.284, 4.284, 4.284, 4.284, 4.721, 5.193, 6.169, 7.129, 7.236, 7.352, 8.397, &
     327  2.163, 2.163, 2.163, 2.163, 2.163, 2.384, 2.622, 3.115, 3.600, 3.654, 3.712, 4.240, &
     328  0.966, 0.966, 0.966, 0.966, 0.966, 1.065, 1.171, 1.392, 1.608, 1.632, 1.658, 1.894 /
     329
     330  DATA cg_MG_6bands/ &
     331     ! bc content = 0.001
     332  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
     333  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
     334  0.643, 0.643, 0.643, 0.654, 0.662, 0.670, 0.677, 0.688, 0.698, 0.704, 0.707, 0.715, &
     335  0.513, 0.513, 0.513, 0.522, 0.530, 0.536, 0.542, 0.552, 0.560, 0.565, 0.568, 0.575, &
     336  0.321, 0.321, 0.321, 0.323, 0.325, 0.327, 0.328, 0.331, 0.333, 0.334, 0.335, 0.337, &
     337  0.153, 0.153, 0.153, 0.149, 0.145, 0.142, 0.139, 0.135, 0.130, 0.128, 0.127, 0.123, &
     338     ! bc content = 0.01
     339  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
     340  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
     341  0.643, 0.643, 0.643, 0.654, 0.662, 0.670, 0.677, 0.688, 0.698, 0.704, 0.707, 0.715, &
     342  0.513, 0.513, 0.513, 0.522, 0.530, 0.536, 0.542, 0.552, 0.560, 0.565, 0.568, 0.575, &
     343  0.321, 0.321, 0.321, 0.323, 0.325, 0.327, 0.328, 0.331, 0.333, 0.334, 0.335, 0.337, &
     344  0.153, 0.153, 0.153, 0.149, 0.145, 0.142, 0.139, 0.135, 0.130, 0.128, 0.127, 0.123, &
     345     ! bc content = 0.02
     346  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
     347  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
     348  0.643, 0.643, 0.643, 0.654, 0.662, 0.670, 0.677, 0.688, 0.698, 0.704, 0.707, 0.715, &
     349  0.513, 0.513, 0.513, 0.522, 0.530, 0.536, 0.542, 0.552, 0.560, 0.565, 0.568, 0.575, &
     350  0.321, 0.321, 0.321, 0.323, 0.325, 0.327, 0.328, 0.331, 0.333, 0.334, 0.335, 0.337, &
     351  0.153, 0.153, 0.153, 0.149, 0.145, 0.142, 0.139, 0.135, 0.130, 0.128, 0.127, 0.123, &
     352     ! bc content = 0.05
     353  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
     354  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
     355  0.643, 0.643, 0.643, 0.654, 0.662, 0.670, 0.677, 0.688, 0.698, 0.704, 0.707, 0.715, &
     356  0.513, 0.513, 0.513, 0.522, 0.530, 0.536, 0.542, 0.552, 0.560, 0.565, 0.568, 0.575, &
     357  0.321, 0.321, 0.321, 0.323, 0.325, 0.327, 0.328, 0.331, 0.333, 0.334, 0.335, 0.337, &
     358  0.153, 0.153, 0.153, 0.149, 0.145, 0.142, 0.139, 0.135, 0.130, 0.128, 0.127, 0.123, &
     359     ! bc content = 0.10
     360  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
     361  0.721, 0.721, 0.721, 0.729, 0.735, 0.741, 0.746, 0.754, 0.762, 0.766, 0.769, 0.775, &
     362  0.643, 0.643, 0.643, 0.654, 0.662, 0.670, 0.677, 0.688, 0.698, 0.704, 0.707, 0.715, &
     363  0.513, 0.513, 0.513, 0.522, 0.530, 0.536, 0.542, 0.552, 0.560, 0.565, 0.568, 0.575, &
     364  0.321, 0.321, 0.321, 0.323, 0.325, 0.327, 0.328, 0.331, 0.333, 0.334, 0.335, 0.337, &
     365  0.153, 0.153, 0.153, 0.149, 0.145, 0.142, 0.139, 0.135, 0.130, 0.128, 0.127, 0.123 /
     366
     367  DATA piz_MG_6bands/ &
     368     ! bc content = 0.001
     369  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
     370  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
     371  0.445, 0.445, 0.445, 0.445, 0.445, 0.521, 0.583, 0.679, 0.741, 0.747, 0.753, 0.798, &
     372  0.394, 0.394, 0.394, 0.394, 0.394, 0.477, 0.545, 0.649, 0.718, 0.724, 0.730, 0.779, &
     373  0.267, 0.267, 0.267, 0.267, 0.267, 0.365, 0.446, 0.571, 0.652, 0.660, 0.667, 0.725, &
     374  0.121, 0.121, 0.121, 0.121, 0.121, 0.139, 0.155, 0.178, 0.193, 0.195, 0.196, 0.207, &
     375     ! bc content = 0.01
     376  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
     377  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
     378  0.445, 0.445, 0.445, 0.445, 0.445, 0.521, 0.583, 0.679, 0.741, 0.747, 0.753, 0.798, &
     379  0.394, 0.394, 0.394, 0.394, 0.394, 0.477, 0.545, 0.649, 0.718, 0.724, 0.730, 0.779, &
     380  0.267, 0.267, 0.267, 0.267, 0.267, 0.365, 0.446, 0.571, 0.652, 0.660, 0.667, 0.725, &
     381  0.121, 0.121, 0.121, 0.121, 0.121, 0.139, 0.155, 0.178, 0.193, 0.195, 0.196, 0.207, &
     382     ! bc content = 0.02
     383  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
     384  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
     385  0.445, 0.445, 0.445, 0.445, 0.445, 0.521, 0.583, 0.679, 0.741, 0.747, 0.753, 0.798, &
     386  0.394, 0.394, 0.394, 0.394, 0.394, 0.477, 0.545, 0.649, 0.718, 0.724, 0.730, 0.779, &
     387  0.267, 0.267, 0.267, 0.267, 0.267, 0.365, 0.446, 0.571, 0.652, 0.660, 0.667, 0.725, &
     388  0.121, 0.121, 0.121, 0.121, 0.121, 0.139, 0.155, 0.178, 0.193, 0.195, 0.196, 0.207, &
     389     ! bc content = 0.05
     390  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
     391  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
     392  0.445, 0.445, 0.445, 0.445, 0.445, 0.521, 0.583, 0.679, 0.741, 0.747, 0.753, 0.798, &
     393  0.394, 0.394, 0.394, 0.394, 0.394, 0.477, 0.545, 0.649, 0.718, 0.724, 0.730, 0.779, &
     394  0.267, 0.267, 0.267, 0.267, 0.267, 0.365, 0.446, 0.571, 0.652, 0.660, 0.667, 0.725, &
     395  0.121, 0.121, 0.121, 0.121, 0.121, 0.139, 0.155, 0.178, 0.193, 0.195, 0.196, 0.207, &
     396     ! bc content = 0.10
     397  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
     398  0.460, 0.460, 0.460, 0.460, 0.460, 0.534, 0.594, 0.688, 0.748, 0.754, 0.760, 0.803, &
     399  0.445, 0.445, 0.445, 0.445, 0.445, 0.521, 0.583, 0.679, 0.741, 0.747, 0.753, 0.798, &
     400  0.394, 0.394, 0.394, 0.394, 0.394, 0.477, 0.545, 0.649, 0.718, 0.724, 0.730, 0.779, &
     401  0.267, 0.267, 0.267, 0.267, 0.267, 0.365, 0.446, 0.571, 0.652, 0.660, 0.667, 0.725, &
     402  0.121, 0.121, 0.121, 0.121, 0.121, 0.139, 0.155, 0.178, 0.193, 0.195, 0.196, 0.207 /
     403
    270404!----BEGINNING OF CALCULATIONS
    271405
    … …  
    342476
    343477  tau_ae(:,:,:,:)=0.
    344   tau_ae_pi(:,:,:,:)=0.
    345478  piz_ae(:,:,:,:)=0.
    346479  cg_ae(:,:,:,:)=0.
     480
     481  tau_ae_pi(:,:,:,:)=0.
     482  piz_ae_pi(:,:,:,:)=0.
     483  cg_ae_pi(:,:,:,:)=0.
    347484   
    348485  DO m=1,nb_aer   ! tau is only computed for each mass
    … …  
    357494        soluble=.TRUE.
    358495        spsol=3
    359         fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD
     496        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for AOD
    360497     ELSEIF  (aerosol_name(m).EQ.id_CSSO4M_phy) THEN
    361498        soluble=.TRUE.
    362499        spsol=4
    363         fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD
     500        fac=1.375    ! (NH4)2-SO4/SO4 132/96 mass conversion factor for AOD
    364501     ELSEIF (aerosol_name(m).EQ.id_SSSSM_phy) THEN
    365502         soluble=.TRUE.
    … …  
    384521     ENDIF
    385522
     523    !--shortname for aerosol index
    386524    id=aerosol_name(m)
    387525
    388526    IF (soluble) THEN
    389527
    390        DO k=1, klev
    391          DO i=1, klon
    392            tmp_var=m_allaer(i,k,spsol)/1.e6*zdh(i,k)*fac
    393            tmp_var_pi=m_allaer_pi(i,k,spsol)/1.e6*zdh(i,k)*fac
    394 
    395            DO inu=1,NSW
    396 
    397              tau_ae2b_int= alpha_aers_6bands(RH_num(i,k),inu,spsol)+ &
    398                            delta(i,k)* (alpha_aers_6bands(RH_num(i,k)+1,inu,spsol) - &
    399                            alpha_aers_6bands(RH_num(i,k),inu,spsol))
     528       !--here we treat the special case of soluble BC internal mixture with Maxwell-Garnett rule
     529       IF (spsol.EQ.1 .AND. flag_bc_internal_mixture) THEN
     530
     531         DO k=1, klev
     532           DO i=1, klon
     533
     534             tmp_var=m_allaer(i,k,spsol)/1.e6*zdh(i,k)*fac
     535             tmp_var_pi=m_allaer_pi(i,k,spsol)/1.e6*zdh(i,k)*fac
     536
     537             ! Calculate the dry BC/(BC+SUL) mass ratio for all (natural+anthropogenic) aerosols
     538             BC_massfra = m_allaer(i,k,id_ASBCM_phy)/(m_allaer(i,k,id_ASBCM_phy)+m_allaer(i,k,id_ASSO4M_phy))
     539
     540             IF (BC_massfra.GE.0.10) THEN
     541               classbc = 5
     542             ELSEIF  (BC_massfra.GE.0.05) THEN
     543               classbc = 4
     544             ELSEIF  (BC_massfra.GE.0.02) THEN
     545               classbc = 3
     546             ELSEIF  (BC_massfra.GE.0.01) THEN
     547               classbc = 2
     548             ELSE
     549               classbc = 1
     550             ENDIF
     551
     552             ! Calculate the dry BC/(BC+SUL) mass ratio for natural aerosols
     553             BC_massfra_pi = m_allaer_pi(i,k,id_ASBCM_phy)/(m_allaer_pi(i,k,id_ASBCM_phy)+m_allaer_pi(i,k,id_ASSO4M_phy))
     554
     555             IF (BC_massfra_pi.GE.0.10) THEN
     556               classbc_pi = 5
     557             ELSEIF  (BC_massfra_pi.GE.0.05) THEN
     558               classbc_pi = 4
     559             ELSEIF  (BC_massfra_pi.GE.0.02) THEN
     560               classbc_pi = 3
     561             ELSEIF  (BC_massfra_pi.GE.0.01) THEN
     562               classbc_pi = 2
     563             ELSE
     564               classbc_pi = 1
     565             ENDIF
     566
     567             DO inu=1,NSW
     568
     569               !--all aerosols
     570               tau_ae2b_int= alpha_MG_6bands(RH_num(i,k),inu,classbc)+                 &
     571                             delta(i,k)* (alpha_MG_6bands(RH_num(i,k)+1,inu,classbc) - &
     572                             alpha_MG_6bands(RH_num(i,k),inu,classbc))
    400573                   
    401              piz_ae2b_int = piz_aers_6bands(RH_num(i,k),inu,spsol) + &
    402                             delta(i,k)* (piz_aers_6bands(RH_num(i,k)+1,inu,spsol) - &
    403                             piz_aers_6bands(RH_num(i,k),inu,spsol))
     574               piz_ae2b_int = piz_MG_6bands(RH_num(i,k),inu,classbc) +                 &
     575                              delta(i,k)* (piz_MG_6bands(RH_num(i,k)+1,inu,classbc) -  &
     576                              piz_MG_6bands(RH_num(i,k),inu,classbc))
    404577                   
    405              cg_ae2b_int = cg_aers_6bands(RH_num(i,k),inu,spsol) + &
    406                            delta(i,k)* (cg_aers_6bands(RH_num(i,k)+1,inu,spsol) - &
    407                            cg_aers_6bands(RH_num(i,k),inu,spsol))
    408 
    409              tau_ae(i,k,id,inu)    = tmp_var*tau_ae2b_int
    410              tau_ae_pi(i,k,id,inu) = tmp_var_pi* tau_ae2b_int
    411              piz_ae(i,k,id,inu)    = piz_ae2b_int
    412              cg_ae(i,k,id,inu)     = cg_ae2b_int
    413                      
     578               cg_ae2b_int = cg_MG_6bands(RH_num(i,k),inu,classbc) +                   &
     579                             delta(i,k)* (cg_MG_6bands(RH_num(i,k)+1,inu,classbc) -    &
     580                             cg_MG_6bands(RH_num(i,k),inu,classbc))
     581
     582               tau_ae(i,k,id,inu)    = tmp_var*tau_ae2b_int
     583               piz_ae(i,k,id,inu)    = piz_ae2b_int
     584               cg_ae(i,k,id,inu)     = cg_ae2b_int
     585
     586               !--natural aerosols
     587               tau_ae2b_int= alpha_MG_6bands(RH_num(i,k),inu,classbc_pi)+                 &
     588                             delta(i,k)* (alpha_MG_6bands(RH_num(i,k)+1,inu,classbc_pi) - &
     589                             alpha_MG_6bands(RH_num(i,k),inu,classbc_pi))
     590                   
     591               piz_ae2b_int = piz_MG_6bands(RH_num(i,k),inu,classbc_pi) +                 &
     592                              delta(i,k)* (piz_MG_6bands(RH_num(i,k)+1,inu,classbc_pi) -  &
     593                              piz_MG_6bands(RH_num(i,k),inu,classbc_pi))
     594                   
     595               cg_ae2b_int = cg_MG_6bands(RH_num(i,k),inu,classbc_pi) +                   &
     596                             delta(i,k)* (cg_MG_6bands(RH_num(i,k)+1,inu,classbc_pi) -    &
     597                             cg_MG_6bands(RH_num(i,k),inu,classbc_pi))
     598
     599               tau_ae_pi(i,k,id,inu) = tmp_var_pi* tau_ae2b_int
     600               piz_ae_pi(i,k,id,inu) = piz_ae2b_int
     601               cg_ae_pi(i,k,id,inu)  = cg_ae2b_int
     602                       
     603             ENDDO
    414604           ENDDO
    415605         ENDDO
    416        ENDDO
     606
     607       !--else treat all other cases of soluble aerosols
     608       ELSE
     609
     610         DO k=1, klev
     611           DO i=1, klon
     612             tmp_var=m_allaer(i,k,spsol)/1.e6*zdh(i,k)*fac
     613             tmp_var_pi=m_allaer_pi(i,k,spsol)/1.e6*zdh(i,k)*fac
     614
     615             DO inu=1,NSW
     616
     617               tau_ae2b_int= alpha_aers_6bands(RH_num(i,k),inu,spsol)+ &
     618                             delta(i,k)* (alpha_aers_6bands(RH_num(i,k)+1,inu,spsol) - &
     619                             alpha_aers_6bands(RH_num(i,k),inu,spsol))
     620                   
     621               piz_ae2b_int = piz_aers_6bands(RH_num(i,k),inu,spsol) + &
     622                            delta(i,k)* (piz_aers_6bands(RH_num(i,k)+1,inu,spsol) - &
     623                              piz_aers_6bands(RH_num(i,k),inu,spsol))
     624                   
     625               cg_ae2b_int = cg_aers_6bands(RH_num(i,k),inu,spsol) + &
     626                             delta(i,k)* (cg_aers_6bands(RH_num(i,k)+1,inu,spsol) - &
     627                             cg_aers_6bands(RH_num(i,k),inu,spsol))
     628
     629               tau_ae(i,k,id,inu)    = tmp_var*tau_ae2b_int
     630               tau_ae_pi(i,k,id,inu) = tmp_var_pi* tau_ae2b_int
     631               piz_ae(i,k,id,inu)    = piz_ae2b_int
     632               cg_ae(i,k,id,inu)     = cg_ae2b_int
     633                       
     634             ENDDO
     635           ENDDO
     636         ENDDO
     637
     638         !--external mixture case for soluble BC
     639         IF (spsol.EQ.1) THEN
     640           piz_ae_pi(:,:,id,:) = piz_ae(:,:,id,:)
     641           cg_ae_pi(:,:,id,:)  = cg_ae(:,:,id,:)
     642         ENDIF
     643
     644       ENDIF
    417645       
    418646     ELSE    ! For all aerosol insoluble components
    … …  
    443671     DO k=1, klev
    444672       DO i=1, klon
    445 !--anthropogenic aerosol
     673!--all (natural + anthropogenic) aerosol
    446674         tau_allaer(i,k,2,inu)=tau_ae(i,k,id_ASSO4M_phy,inu)+tau_ae(i,k,id_CSSO4M_phy,inu)+ &
    447675                               tau_ae(i,k,id_ASBCM_phy,inu)+tau_ae(i,k,id_AIBCM_phy,inu)+   &
    … …  
    449677                               tau_ae(i,k,id_ASSSM_phy,inu)+tau_ae(i,k,id_CSSSM_phy,inu)+   &
    450678                               tau_ae(i,k,id_SSSSM_phy,inu)+ tau_ae(i,k,id_CIDUSTM_phy,inu)
    451          tau_allaer(i,k,2,inu)=MAX(tau_allaer(i,k,2,inu),1e-15)
    452 
    453          piz_allaer(i,k,2,inu)=(tau_ae(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)+ &
    454                                 tau_ae(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)+ &
    455                                 tau_ae(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)+ &
    456                                 tau_ae(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)+ &
    457                                 tau_ae(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)+ &
    458                                 tau_ae(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)+ &
    459                                 tau_ae(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)+ &
    460                                 tau_ae(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)+ &
    461                                 tau_ae(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)+ &
     679         tau_allaer(i,k,2,inu)=MAX(tau_allaer(i,k,2,inu),tau_min)
     680
     681         piz_allaer(i,k,2,inu)=(tau_ae(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)+   &
     682                                tau_ae(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)+   &
     683                                tau_ae(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)+     &
     684                                tau_ae(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)+     &
     685                                tau_ae(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)+   &
     686                                tau_ae(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)+   &
     687                                tau_ae(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)+     &
     688                                tau_ae(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)+     &
     689                                tau_ae(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)+     &
    462690                                tau_ae(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)) &
    463691                                /tau_allaer(i,k,2,inu)
    464          piz_allaer(i,k,2,inu)=MAX(piz_allaer(i,k,2,inu),0.01)
     692         piz_allaer(i,k,2,inu)=MIN(MAX(piz_allaer(i,k,2,inu),0.01),1.0)
     693         IF (tau_allaer(i,k,2,inu).LE.tau_min) piz_allaer(i,k,2,inu)=1.0
    465694
    466695         cg_allaer(i,k,2,inu)=(tau_ae(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)*cg_ae(i,k,id_ASSO4M_phy,inu)+ &
    467696                               tau_ae(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)*cg_ae(i,k,id_CSSO4M_phy,inu)+ &
    468                                tau_ae(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)*cg_ae(i,k,id_ASBCM_phy,inu)+ &
    469                                tau_ae(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)*cg_ae(i,k,id_AIBCM_phy,inu)+ &
     697                               tau_ae(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)*cg_ae(i,k,id_ASBCM_phy,inu)+    &
     698                               tau_ae(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)*cg_ae(i,k,id_AIBCM_phy,inu)+    &
    470699                               tau_ae(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)*cg_ae(i,k,id_ASPOMM_phy,inu)+ &
    471700                               tau_ae(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)*cg_ae(i,k,id_AIPOMM_phy,inu)+ &
    472                                tau_ae(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)*cg_ae(i,k,id_ASSSM_phy,inu)+ &
    473                                tau_ae(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)*cg_ae(i,k,id_CSSSM_phy,inu)+ &
    474                                tau_ae(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)*cg_ae(i,k,id_SSSSM_phy,inu)+ &
     701                               tau_ae(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)*cg_ae(i,k,id_ASSSM_phy,inu)+    &
     702                               tau_ae(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)*cg_ae(i,k,id_CSSSM_phy,inu)+    &
     703                               tau_ae(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)*cg_ae(i,k,id_SSSSM_phy,inu)+    &
    475704                               tau_ae(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)*cg_ae(i,k,id_CIDUSTM_phy,inu))/ &
    476705                               (tau_allaer(i,k,2,inu)*piz_allaer(i,k,2,inu))
     706         cg_allaer(i,k,2,inu)=MIN(MAX(cg_allaer(i,k,2,inu),0.0),1.0)
    477707
    478708!--natural aerosol
     709!--ASBCM aerosols take _pi value because of internal mixture option
    479710         tau_allaer(i,k,1,inu)=tau_ae_pi(i,k,id_ASSO4M_phy,inu)+tau_ae_pi(i,k,id_CSSO4M_phy,inu)+ &
    480711                               tau_ae_pi(i,k,id_ASBCM_phy,inu)+tau_ae_pi(i,k,id_AIBCM_phy,inu)+   &
    … …  
    482713                               tau_ae_pi(i,k,id_ASSSM_phy,inu)+tau_ae_pi(i,k,id_CSSSM_phy,inu)+   &
    483714                               tau_ae_pi(i,k,id_SSSSM_phy,inu)+ tau_ae_pi(i,k,id_CIDUSTM_phy,inu)
    484          tau_allaer(i,k,1,inu)=MAX(tau_allaer(i,k,1,inu),1e-15)
    485 
    486          piz_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)+ &
    487                                 tau_ae_pi(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)+ &
    488                                 tau_ae_pi(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)+ &
    489                                 tau_ae_pi(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)+ &
    490                                 tau_ae_pi(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)+ &
    491                                 tau_ae_pi(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)+ &
    492                                 tau_ae_pi(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)+ &
    493                                 tau_ae_pi(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)+ &
    494                                 tau_ae_pi(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)+ &
     715         tau_allaer(i,k,1,inu)=MAX(tau_allaer(i,k,1,inu),tau_min)
     716
     717         piz_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)+   &
     718                                tau_ae_pi(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)+   &
     719                                tau_ae_pi(i,k,id_ASBCM_phy,inu)*piz_ae_pi(i,k,id_ASBCM_phy,inu)+ &
     720                                tau_ae_pi(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)+     &
     721                                tau_ae_pi(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)+   &
     722                                tau_ae_pi(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)+   &
     723                                tau_ae_pi(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)+     &
     724                                tau_ae_pi(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)+     &
     725                                tau_ae_pi(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)+     &
    495726                                tau_ae_pi(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)) &
    496727                                /tau_allaer(i,k,1,inu)
    497          piz_allaer(i,k,1,inu)=MAX(piz_allaer(i,k,1,inu),0.01)
    498 
    499          cg_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)*cg_ae(i,k,id_ASSO4M_phy,inu)+ &
    500                                tau_ae_pi(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)*cg_ae(i,k,id_CSSO4M_phy,inu)+ &
    501                                tau_ae_pi(i,k,id_ASBCM_phy,inu)*piz_ae(i,k,id_ASBCM_phy,inu)*cg_ae(i,k,id_ASBCM_phy,inu)+ &
    502                                tau_ae_pi(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)*cg_ae(i,k,id_AIBCM_phy,inu)+ &
    503                                tau_ae_pi(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)*cg_ae(i,k,id_ASPOMM_phy,inu)+ &
    504                                tau_ae_pi(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)*cg_ae(i,k,id_AIPOMM_phy,inu)+ &
    505                                tau_ae_pi(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)*cg_ae(i,k,id_ASSSM_phy,inu)+ &
    506                                tau_ae_pi(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)*cg_ae(i,k,id_CSSSM_phy,inu)+ &
    507                                tau_ae_pi(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)*cg_ae(i,k,id_SSSSM_phy,inu)+ &
     728         piz_allaer(i,k,1,inu)=MIN(MAX(piz_allaer(i,k,1,inu),0.01),1.0)
     729         IF (tau_allaer(i,k,1,inu).LE.tau_min) piz_allaer(i,k,1,inu)=1.0
     730
     731         cg_allaer(i,k,1,inu)=(tau_ae_pi(i,k,id_ASSO4M_phy,inu)*piz_ae(i,k,id_ASSO4M_phy,inu)*cg_ae(i,k,id_ASSO4M_phy,inu)+    &
     732                               tau_ae_pi(i,k,id_CSSO4M_phy,inu)*piz_ae(i,k,id_CSSO4M_phy,inu)*cg_ae(i,k,id_CSSO4M_phy,inu)+    &
     733                               tau_ae_pi(i,k,id_ASBCM_phy,inu)*piz_ae_pi(i,k,id_ASBCM_phy,inu)*cg_ae_pi(i,k,id_ASBCM_phy,inu)+ &
     734                               tau_ae_pi(i,k,id_AIBCM_phy,inu)*piz_ae(i,k,id_AIBCM_phy,inu)*cg_ae(i,k,id_AIBCM_phy,inu)+       &
     735                               tau_ae_pi(i,k,id_ASPOMM_phy,inu)*piz_ae(i,k,id_ASPOMM_phy,inu)*cg_ae(i,k,id_ASPOMM_phy,inu)+    &
     736                               tau_ae_pi(i,k,id_AIPOMM_phy,inu)*piz_ae(i,k,id_AIPOMM_phy,inu)*cg_ae(i,k,id_AIPOMM_phy,inu)+    &
     737                               tau_ae_pi(i,k,id_ASSSM_phy,inu)*piz_ae(i,k,id_ASSSM_phy,inu)*cg_ae(i,k,id_ASSSM_phy,inu)+       &
     738                               tau_ae_pi(i,k,id_CSSSM_phy,inu)*piz_ae(i,k,id_CSSSM_phy,inu)*cg_ae(i,k,id_CSSSM_phy,inu)+       &
     739                               tau_ae_pi(i,k,id_SSSSM_phy,inu)*piz_ae(i,k,id_SSSSM_phy,inu)*cg_ae(i,k,id_SSSSM_phy,inu)+       &
    508740                               tau_ae_pi(i,k,id_CIDUSTM_phy,inu)*piz_ae(i,k,id_CIDUSTM_phy,inu)*cg_ae(i,k,id_CIDUSTM_phy,inu))/ &
    509741                               (tau_allaer(i,k,1,inu)*piz_allaer(i,k,1,inu))
     742         cg_allaer(i,k,1,inu)=MIN(MAX(cg_allaer(i,k,1,inu),0.0),1.0)
    510743
    511744        ENDDO
    … …  
    513746    ENDDO
    514747   
    515 !--waveband 2 and all aerosol
     748!--waveband 2 and all aerosol (third index = 2)
    516749  inu=2
    517750  DO i=1, klon

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

    r2542 r2669  
    22!
    33SUBROUTINE readaerosol_optic_rrtm(debut, aerosol_couple,  &
    4      new_aod, flag_aerosol, itap, rjourvrai, &
     4     new_aod, flag_aerosol, flag_bc_internal_mixture, itap, rjourvrai, &
    55     pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
    66     tr_seri, mass_solu_aero, mass_solu_aero_pi, &
    … …  
    3232  LOGICAL, INTENT(IN)                      :: new_aod
    3333  INTEGER, INTENT(IN)                      :: flag_aerosol
     34  LOGICAL, INTENT(IN)                      :: flag_bc_internal_mixture
    3435  INTEGER, INTENT(IN)                      :: itap
    3536  REAL, INTENT(IN)                         :: rjourvrai
    … …  
    303304!--new aerosol properties
    304305  ! aeropt_6bands for rrtm
    305   CALL aeropt_6bands_rrtm( &
    306        pdel, m_allaer, pdtphys, rhcl, &
    307        tau_aero, piz_aero, cg_aero,   &
    308        m_allaer_pi, flag_aerosol, &
    309        zrho )
     306  CALL aeropt_6bands_rrtm(          &
     307       pdel, m_allaer, rhcl,        &
     308       tau_aero, piz_aero, cg_aero, &
     309       m_allaer_pi, flag_aerosol,   &
     310       flag_bc_internal_mixture, zrho )
    310311
    311312  ! aeropt_5wv only for validation and diagnostics
    312   CALL aeropt_5wv_rrtm(                    &
    313        pdel, m_allaer,                &
    314        pdtphys, rhcl, aerindex,       &
    315        flag_aerosol, pplay, t_seri,   &
     313  CALL aeropt_5wv_rrtm(              &
     314       pdel, m_allaer,               &
     315       rhcl, aerindex,               &
     316       flag_aerosol,                 &
     317       flag_bc_internal_mixture,     &
     318       pplay, t_seri,                &
    316319       tausum_aero, tau3d_aero )
    317320
Note: See TracChangeset for help on using the changeset viewer.