oi_cacsts.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!option! -O nomove
!****---------------------------------------------------------------------------
!****   CACSTS : INITIALIZES THE SURFACE FIELDS
!****   ------
!****  Auteurs :   CB 01/91, BU 05/92, VC 05/93, DG 03/94, PA 09/95, DG 05/96
!****      Mod : E. Bazile 01/97 Subtraction of the mean temperature and/or
!****                            humidity bias for the increments used for the 
!****                            analysis of the soil water
!****      Mod : D. Giard  03/99 ACSOL -> ACSOLW
!****            E. Bazile , F. Bouyssel : logical LLLACW is replaced by a 
!****                       continuous function ZDACW (LSOLV).
!****      Mod : F. Taillefer 09/02 : update of surface constants according to SST
!****      Mod : F. Bouyssel 02/04 : threshold using the solar zenithal angle
!****      Mod : E. Bazile 01/2007 : Parameter for the correction PSNS and WPI
!        M.Hamrud      01-Jul-2006  Revised surface fields
!        A.Trojakova   27-Jun-2007 bugfixing ZV10M (surface pointers)
!        F. Bouyssel    27-Mar-2011  Use of REPS2 instead of REPS3 for ZNEI
!****---------------------------------------------------------------------------
!
SUBROUTINE oi_cacsts(KNBPT,PT2INC,PH2INC,PWGINC,PWS_O,                      &
                     kdat,ksssss,                                           &
                     ptp,pwp,ptl,psns,pts,pws,                              &
                     ptcls,phcls,pucls,pvcls,psstc,pwpinc1,pwpinc2,pwpinc3, &
                     pt2mbias,ph2mbias,                                     &
                     prrcl,prrsl,prrcn,prrsn,patmneb,pevap,pevaptr,         &
                     pitm,pveg,palbf,pemisf,pz0f,                           &
                     piveg,parg,pd2,psab,plai,prsmin,pz0h,                  &
                     ptsc,ptpc,pwsc,pwpc,psnc,pgelat,pgelam,pgemu)  
!
!****---------------------------------------------------------------------------
!**  AIM : INITIALIZES THE PRONOSTIC SURFACE FIELDS
!**  ---
!**  SEQUENCE D'APPEL :
!**  ----------------
!**        CALL CACSTS(....)

!**  INPUT ARGUMENTS :
!**  ------------------        
!**                               
!**        - EXPLICIT - 
!**                      KNBPT  : real number of treated points
!**                      PT2INC : analysis increment of T2m
!**                      PH2INC : analysis increment of Hu2m
!**                      PSP_SB,PSP_SG,PSP_RR,PSD_VF,PSD_VV,PSD_VX,PSP_CI,PSP_X2    : 
!**                      buffer of pgd analysis fields 
!**                      PGELAM, PGELAT, PGEMU : geographical coordinates

!**  OUTPUT ARGUMENTS : 
!**  -------------------
!**  EXTERN : CAVEGI (FCTVEG) - ACSOLW - TSL
!**  --------

!**  ALGORITHM :  - INITIALIZES THE SURFACE TEMPERATURE.
!**  ----------   - INITIALIZES THE DEEP TEMPERATURE.
!**               - INITIALIZES THE SURFACE WATER TANK.
!**               - INITIALIZES THE DEEP WATER TANK.
!**               - CORRECTS THE SNOW AMOUNT.
!***-----------------------------------------------------------------
!
USE modd_csts,       ONLY : xg, xtt, xrholw, xday
USE modd_surf_par,   ONLY : xundef
USE modd_assim
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_oi_cavegi
USE modi_oi_acsolw
USE modi_oi_jacobians
USE modi_oi_tsl
USE modi_oi_fctveg
USE modi_oi_kalman_gain
!
IMPLICIT NONE
!
INTEGER,INTENT(IN)    :: knbpt, kdat, ksssss
!
REAL   ,INTENT(IN)    :: pt2inc(knbpt) 
REAL   ,INTENT(IN)    :: ph2inc(knbpt) 
REAL   ,INTENT(IN)    :: pwginc(knbpt)
REAL   ,INTENT(IN)    :: pws_o(knbpt)
REAL   ,INTENT(INOUT) :: ptp(knbpt)
REAL   ,INTENT(INOUT) :: pwp(knbpt)
REAL   ,INTENT(INOUT) :: ptl(knbpt)
REAL   ,INTENT(INOUT) :: psns(knbpt)
REAL   ,INTENT(INOUT) :: pts(knbpt)
REAL   ,INTENT(INOUT) :: pws(knbpt)
REAL   ,INTENT(INOUT) :: ptcls(knbpt) 
REAL   ,INTENT(INOUT) :: phcls(knbpt) 
REAL   ,INTENT(INOUT) :: pucls(knbpt)
REAL   ,INTENT(INOUT) :: pvcls(knbpt)
REAL   ,INTENT(INOUT) :: psstc(knbpt)
REAL   ,INTENT(INOUT) :: pwpinc1(knbpt)
REAL   ,INTENT(INOUT) :: pwpinc2(knbpt)
REAL   ,INTENT(INOUT) :: pwpinc3(knbpt)
REAL   ,INTENT(INOUT) :: pt2mbias(knbpt)
REAL   ,INTENT(INOUT) :: ph2mbias(knbpt)
REAL   ,INTENT(IN)    :: prrcl(knbpt)
REAL   ,INTENT(IN)    :: prrsl(knbpt)
REAL   ,INTENT(IN)    :: prrcn(knbpt)
REAL   ,INTENT(IN)    :: prrsn(knbpt)
REAL   ,INTENT(IN)    :: patmneb(knbpt)
REAL   ,INTENT(IN)    :: pevap(knbpt)
REAL   ,INTENT(IN)    :: pevaptr(knbpt)
REAL   ,INTENT(IN)    :: pitm(knbpt) 
REAL   ,INTENT(IN) :: pveg(knbpt) 
REAL   ,INTENT(INOUT) :: palbf(knbpt)
REAL   ,INTENT(INOUT) :: pemisf(knbpt)
REAL   ,INTENT(INOUT) :: pz0f(knbpt)
REAL   ,INTENT(INOUT) :: piveg(knbpt)
REAL   ,INTENT(INOUT) :: parg(knbpt)
REAL   ,INTENT(INOUT) :: pd2(knbpt)
REAL   ,INTENT(INOUT) :: psab(knbpt) 
REAL   ,INTENT(INOUT) :: plai(knbpt)
REAL   ,INTENT(INOUT) :: prsmin(knbpt)
REAL   ,INTENT(INOUT) :: pz0h(knbpt)
REAL   ,INTENT(IN)    :: ptsc(knbpt)
REAL   ,INTENT(IN)    :: ptpc(knbpt)
REAL   ,INTENT(IN)    :: pwsc(knbpt)
REAL   ,INTENT(IN)    :: pwpc(knbpt)
REAL   ,INTENT(IN)    :: psnc(knbpt) 
REAL   ,INTENT(IN)    :: pgelat(knbpt) 
REAL   ,INTENT(IN)    :: pgelam(knbpt) 
REAL   ,INTENT(IN)    :: pgemu(knbpt)
!
REAL, DIMENSION(24) :: zvgat1,zvgat2,zvgat3,zvgbt1,zvgbt2,zvgbt3,zvgct1,zvgct2
REAL, DIMENSION(24) :: zvgah1,zvgah2,zvgah3,zvgbh1,zvgbh2,zvgbh3,zvgch1,zvgch2
REAL, DIMENSION(24) :: zsigt2mp,zsighp2
!
REAL, DIMENSION(KNBPT) :: ziveg
REAL, DIMENSION(KNBPT) :: zwfc, zwpmx, zwsat, zwsmx, zwwilt
REAL, DIMENSION(KNBPT) :: zdwg_dwg, zdwg_dw2
!
REAL :: zechgu, znei, zcli, zpd, zclimca
REAL :: ztsc, ztpc, zwsc, zwpc, zsnc
REAL :: zv10m, zprecip, zwpi, zdacw, zdacw2, zmu0, zmu0m
!
REAL :: zvgst,zvgsh,zvgpt1,zvgph1,zvgpt2,zvgph2,zg1,zg2,zg3,zg4
!
REAL :: zzt, zzh, zlaisrs, zteff, zheff
REAL :: zcwph, zcwpt, zt2d, zh2d
REAL :: zwsd, zwpd, zwpdx
!
REAL :: zwsa, zwsmin, zwpa, zwpmin
REAL :: zgel, zsna, zmsn, zk1, zk2
!
INTEGER :: ih, jrof
LOGICAL :: gsgobs
!
REAL(KIND=JPRB) :: zhook_handle
!
!--------------------------------------------------------------------
IF (lhook) CALL dr_hook('OI_CACSTS',0,zhook_handle)
!
zechgu = REAL(NECHGU) * 3600.
!
!**  1.1 Initialization of raw polynomials and reference fields.   
!
 CALL oi_cavegi(zvgat1,zvgat2,zvgat3,zvgbt1,zvgbt2,zvgbt3,zvgct1,zvgct2, &
                zvgah1,zvgah2,zvgah3,zvgbh1,zvgbh2,zvgbh3,zvgch1,zvgch2, &
                zsigt2mp,zsighp2,gsgobs) 
!
!
!*   1.2  Initialization of intermediate variables
!
DO jrof = 1,knbpt
  ziveg(jrof) = anint(piveg(jrof))
ENDDO
!
 CALL oi_acsolw(1, knbpt,                        &
                 parg, pd2, pws, ziveg, psab,     &
                 lldhmt,                          &
                 zwfc, zwpmx, zwsat, zwsmx, zwwilt)  
!
!* 1.3. Analytical Jacobians for WG assimilation
!
 CALL oi_jacobians(knbpt,pws_o,psab,parg,pd2,pwp,zdwg_dwg,zdwg_dw2) 
!
!**---------------------------------------------------------------------
!**  - 2 - Calculation of analysed fields.

DO jrof = 1,knbpt

! surface analysis or  de surface ou call back to climatology, on earth
!    IF (PITM(JROF) > 0.5.AND.(RCLIMCA >= 0.0).AND.PWS(JROF)/=XUNDEF) THEN
  IF ( pws(jrof)/=xundef ) THEN
!
! storage of forecast fields
    znei = max(0.0,psns(jrof)/(psns(jrof)+xwcrin))
! update of climatological fields    
    zcli = xrclimca /(1.0+xrclimn*znei)

    IF ( .NOT. lclim ) THEN
      ztsc = pts(jrof)
      ztpc = ptp(jrof)
      zwsc = pws(jrof)
      zwpc = pwp(jrof)
      zsnc = psns(jrof)
    ELSE
      ztsc = ptsc(jrof)
      ztpc = ptpc(jrof)
      zwsc = pwsc(jrof) * zwsmx(jrof)
      zwpc = pwpc(jrof) * zwpmx(jrof)
      zsnc = psnc(jrof)
    ENDIF
    
!-----------------------------------------------------------------------------------
!
!*   2.1  Temperature analysis
!
! transfer of 2m temperature increment on Ts and Tp with damping
!
    IF ( nneigt<=0 .OR. znei<xreps2 ) THEN
      zpd = 1.0
    ELSEIF ( xsneigt<xreps3 ) THEN
      zpd = 0.0
    ELSE
      zpd = (1.0-min(znei,xsneigt)/xsneigt)**nneigt
    ENDIF

    pts(jrof) =  pts(jrof)  + pt2inc(jrof)*zpd
    ptp(jrof) =  ptp(jrof)  + pt2inc(jrof)*zpd/(xsodelx(1)/xsodelx(0))

! Call back of Ts
    zclimca = xrclimts * zcli
    pts(jrof) = (1.0-zclimca)*pts(jrof) + zclimca*ztsc

! Call back of Tp
    zclimca = xrclimtp * xrclimca
    ptp(jrof) = (1.0-zclimca)*ptp(jrof )+ zclimca*ztpc

!-----------------------------------------------------------------------------------
!
!*   2.2  Initializations for the surface analysis
!
! local conditions for the effective analysis of surface fields
! calculation of the useful local solar time
!
    CALL oi_tsl(kdat,ksssss,pgelat(jrof),pgelam(jrof),zmu0,zmu0m,ih)
!
    zv10m = sqrt(pucls(jrof)**2+pvcls(jrof)**2)
!
    zprecip = max(0.,prrcl(jrof))+ max(0.,prrsl(jrof)) &
            + max(0.,prrcn(jrof))+ max(0.,prrsn(jrof))  
!
    IF (lfgel) THEN
      zwpi = ptl(jrof)
    ELSE
      zwpi = 0.0
    ENDIF
!
! Surface water forcing to the superficial reservoir 
!
    zdacw =  min(1.0,max(0.0,abs(REAL(nint(ziveg(jrof))-ntvgla))) )  &
           * min(1.0,max(0.0,REAL(ih)))                              &
           * min(1.0,max(0.0,REAL(nidj)/REAL(nmindj)))               &
           * min(1.0,max(0.0,1.0-zv10m/(xv10mx+xreps3)))             &
           * min(1.0,max(0.0,1.0-zprecip/(xsprecip+xreps3)))         &
           * min(1.0,max(0.0,1.0-zwpi/xsice))  
!
! coefficients : depend on the solar zenithal angle
!
    IF ( xsmu0>xreps3 ) THEN
      zpd = 0.5 * (1.0+tanh(xsmu0*(zmu0m-0.5)))
    ELSE
      zpd = 1.0
    ENDIF
    zdacw = zdacw * zpd
!
! coefficients : depend on the surface evaporation
!
!*    Threshold of min. evaporation for W analysis (SEVAP en mm/day)
    IF ( xsevap>xreps3 ) THEN
      zpd = min(1.0,max(0.0,pevap(jrof)/(-xsevap/xday)))
    ELSE
      zpd = 1.0
    ENDIF
    zdacw = zdacw * zpd
!
! coefficients : depend on the nebulosity
!
    IF ( xanebul>xreps3 ) THEN
      zpd = 1.0 - xanebul*(patmneb(jrof)/zechgu)**nnebul
    ELSE
      zpd = 1.0
    ENDIF
    zdacw = zdacw * zpd
!
! coefficients : depend on the snow cover

    IF ( nneigw<=0 .OR. znei<xreps2 ) THEN
      zpd = 1.0
    ELSEIF ( xsneigw<xreps3 ) THEN
      zpd = 0.0
    ELSE
      zpd = ( 1.0 - min(znei,xsneigw)/xsneigw)**nneigw
    ENDIF
    zdacw = zdacw * zpd
!
    zdacw2 = min(1.0,max(0.0,1.0-(zprecip+abs(pevap(jrof)))/(xsprecip2+xreps3)))  &
           * min(1.0,max(0.0,1.0-zwpi/xsice))
!
    zdacw2 = zdacw2 * zpd


!*   2.3  Humidity analysis by optimal interpolation for ISBA


! coefficients : mainly depend on vegetation
!
!  fctveg.h 
!****---------------Calculation of ZWSD and ZWPD------------------------------------
!
    CALL oi_fctveg(ih,pveg(jrof),                                            &
                   zvgat1,zvgat2,zvgat3,zvgbt1,zvgbt2,zvgbt3,zvgct1,zvgct2,  &
                   zvgah1,zvgah2,zvgah3,zvgbh1,zvgbh2,zvgbh3,zvgch1,zvgch2,  &
                   zsigt2mp,zsighp2,                                         &
                   zg1,zg2,zg3,zg4,                                          &
                   zvgst,zvgsh,zvgpt1,zvgph1,zvgpt2,zvgph2)  
!
! coefficients : depend on the observation errors
! nb - in our case GSGOBS=.F.
!
    IF ( gsgobs ) THEN
      zzt = zg1 / zg2
      zzh = zg3 / zg4 
    ELSE
      zzt = 1.0
      zzh = 1.0
    ENDIF
    
! coefficients : depend on the texture
!
    zpd = (zwfc(jrof)-zwwilt(jrof))/xadwr

! calculation of raw increments for ws=Ws/ds/ro, wp=Wp/dp/ro
! final coefficients
!
    zzt = zzt * zpd * zdacw
    zzh = zzh * zpd * zdacw
!
    zvgst = zvgst * zzt
    zvgsh = zvgsh * zzh

    zlaisrs = plai(jrof)/max(1.0,prsmin(jrof))    
    zcwpt   = ( zvgpt1 + zlaisrs*zvgpt2 ) * zzt
    zcwph   = ( zvgph1 + zlaisrs*zvgph2 ) * zzh
!
! possible limitation of increments for T2m and H2m
! limitation of the absolute value of the increments
!
    zt2d = pt2inc(jrof)
    zh2d = ph2inc(jrof)
    IF (xsigt2mo < 0.0) zt2d=max(xsigt2mo,min(-xsigt2mo,zt2d))
    IF (xsigh2mo < 0.0) zh2d=max(xsigh2mo,min(-xsigh2mo,zh2d))

! removal of the mean bias
! subtraction of mean biais if SCOEF(T/H) <> 1
!
    pt2mbias(jrof) = pt2mbias(jrof)*(1.0-xscoeft) + zt2d*xscoeft
    ph2mbias(jrof) = ph2mbias(jrof)*(1.0-xscoefh) + zh2d*xscoefh

! if the current bias is lower than the mean bias, it's set to zero
!                IF (ABS(ZT2D).LT.ABS(PSP_CI(JROF,YSP_CI%YCI(12)%MP0)) ZTEFF = 0.
!                IF (ABS(ZH2D).LT.ABS(PSP_CI(JROF,YSP_CI%YCI(13)%MP0)) ZHEFF = 0.
! if the current bias is lower than the effective bias, it's kept
!
    IF ( xscoeft/= 0.0 .OR. xscoefh/=0.0 ) THEN
      zteff = zt2d - pt2mbias(jrof)
      zheff = zh2d - ph2mbias(jrof)
      IF (abs(zt2d) < abs(zteff)) zteff = zt2d
      IF (abs(zh2d) < abs(zheff)) zheff = zh2d
      zt2d = zteff
      zh2d = zheff
    ENDIF

! raw increments
!
    IF ( lobs2m .AND. (.NOT.lobswg .OR. pwginc(jrof)==0.0) ) THEN
      zwsd = xrscaldw * (zvgst*zt2d + zvgsh*zh2d)
      zwpd = xrscaldw * (zcwpt*zt2d + zcwph*zh2d)
    ELSEIF ( lobswg ) THEN
      CALL oi_kalman_gain(zdwg_dwg(jrof),zdwg_dw2(jrof),pd2(jrof),zk1,zk2)
      zwsd = zk1*zdacw2*pwginc(jrof)
      zwpd = zk2*zdacw2*pwginc(jrof)
    ELSE
      zwsd = 0.0
      zwpd = 0.0
    ENDIF

! limitations on the corrections
! no ws analysis if no evaporation on bare soil
!
    IF (pevap(jrof)-pevaptr(jrof)>= 0.0 .AND. .NOT.lobswg) zwsd = 0.0

!===============================================================
! Lower limit for Wp set to  Wwilt instead of veg*Wwilt
!===============================================================
!        ZZVEG = 1.0
!=========================WS et WP================================

! analysis of wp limited to assure veg*wwilt <= wp <= SWFC*wfc
    IF ( limveg ) CALL get_zw(zwpd,pwp(jrof),pd2(jrof))

    ! smoothing of analysis increments of wp
    IF ( lissew ) THEN
      zwpdx = zwpd
      IF ( nlissew >= 3 ) THEN
        zwpd = 0.25*(pwpinc3(jrof)+pwpinc2(jrof)+pwpinc1(jrof)+zwpdx)    
      ELSE
        zwpd = 0.0
      ENDIF
      IF ( nlissew >= 2 ) pwpinc3(jrof) = pwpinc2(jrof)
      IF ( nlissew >= 1 ) pwpinc2(jrof) = pwpinc1(jrof)
      pwpinc1(jrof) = zwpdx
    ENDIF

! analysis of ws limited to assure veg*wwilt <= ws <= SWFC*wfc
    IF ( limveg ) CALL get_zw(zwsd,pws(jrof),xrd1)

! transfer of the increments on Ws, Wp
!
    zwsa = pws(jrof) + zwsd * xrd1 * xrholw
    zwsmin = xreps1 * xrd1 * xrholw
    pws(jrof) = max(zwsmin,min(zwsmx(jrof),zwsa)) 

! minimal total water contents
!
    zwpa = pwp(jrof) + zwpd * pd2(jrof) * xrholw
    zwpmin = max(pws(jrof), xreps1 * pd2(jrof) * xrholw)
    pwp(jrof) = max(zwpmin,min(zwpmx(jrof),zwpa))

! Call back of Ws
    zclimca = xrclimws * zcli
    zclimca = zclimca*pveg(jrof) + min(1.0,xrclimv*zclimca)*(1.0-pveg(jrof))
    pws(jrof) = (1.0-zclimca)*pws(jrof) + zclimca*zwsc

! Call back of Wp
    zclimca = xrclimwp * zcli
    zclimca = zclimca*pveg(jrof) + min(1.0,xrclimv*zclimca)*(1.0-pveg(jrof))
    IF ( lfgel ) THEN
      zgel = zwpi / max(pwp(jrof)+zwpi,xreps3)
      zwpc = zwpc * (1.0 - max(0.0,min(1.0,zgel)))
      zwpc = max(zwpmin,zwpc)
    ENDIF
    pwp(jrof) = (1.0-zclimca)*pwp(jrof) + zclimca*zwpc

! call back of Sn with a possible correction for melting
!
    zsna  = (1.0-xrclimca)*psns(jrof) + xrclimca*zsnc
    zmsn  = max(0.0, xrsnsa/21600.*zechgu*(ptcls(jrof)-xtt))**xrsnsb
    psns(jrof) = max(zsna-zmsn,0.0)

    IF (lfgel) THEN
      zmsn = max(0.0, xrwpia/21600. * zechgu*(ptcls(jrof)-xtt))**xrwpib
      ptl(jrof)=max(zwpi-zmsn ,0.0)
      pwp(jrof)=pwp(jrof)-ptl(jrof)+zwpi
    ENDIF

!*   2.5  Call back of SST, on sea

!  ELSEIF ( PITM(JROF) <= 0.5 .AND. RCLISST /= 0. .AND. LCLIM ) THEN  
!    PTS(JROF) = (1.0-RCLISST)*PTS(JROF) + RCLISST *PSSTC(JROF)    
!    PTP(JROF) = PTS(JROF)
!   PWS(JROF) = XUNDEF
!   PWP(JROF) = XUNDEF
!   PTL(JROF) = 0.0
  ENDIF

!*   2.6  Update of surface constants on sea, functions of ice field

  IF ( pitm(jrof) <= 0.5 ) THEN
    IF ( pts(jrof) <= xtmergl ) THEN
      palbf(jrof) = xsalbb
      pemisf(jrof) = xsemib
      pz0f(jrof) = xszz0b*xg
      pz0h(jrof) = xrzhz0g * xszz0b*xg
    ELSE
      palbf(jrof) = xsalbm
      pemisf(jrof) = xsemim
    ENDIF
  ENDIF

ENDDO

IF (lhook) CALL dr_hook('OI_CACSTS',1,zhook_handle)

 CONTAINS

SUBROUTINE get_zw(PWD,PIW,PD)
!
REAL, INTENT(INOUT) :: pwd
REAL, INTENT(IN) :: piw, pd
REAL :: zwr, zwd1, zwd2
!
zwr = piw / (pd*xrholw)
!
IF ( zwr > zwfc(jrof)*xswfc ) THEN
  IF ( lhumid ) THEN
    pwd = min(0.0,pwd)
  ELSE
    pwd = 0.0
  ENDIF
ELSEIF ( zwr < zwwilt(jrof)*pveg(jrof)) THEN
  IF (lhumid) THEN
    pwd = max(0.0,pwd)
  ELSE
    pwd = 0.0
  ENDIF
ELSE
  zwd1 = zwwilt(jrof)*pveg(jrof) - zwr
  zwd2 = zwfc(jrof)*xswfc   - zwr
  pwd = max(zwd1,min(zwd2,pwd))
ENDIF

END SUBROUTINE
!
!**---------------------------------------------------------------------
END SUBROUTINE oi_cacsts