ol_time_interp_atm.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.
!     ######spl
SUBROUTINE ol_time_interp_atm (KSURF_STEP,KNB_ATM,KSIZE_OMP,             &
                               pta1,pta2,pqa1,pqa2,pwind1,pwind2,        &
                               pdir_sw1,pdir_sw2,psca_sw1,psca_sw2,      &
                               plw1,plw2,psnow2,prain2,                  &
                               pps1,pps2,pco21,pco22,pdir1,pdir2         )  
!**************************************************************************
!
!!    PURPOSE
!!    -------
!        Time interpolation of the atmospheric forcing
!        So far, it is a simple linear interpolation.
!        More complex interpolation may be added, especially for the atmospheric
!        radiation (Option to use).
!        Output are in the module 
!!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!
!!    AUTHOR
!!    ------
!!      F. Habets   *Meteo France*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    06/2003
!
USE modd_csts,       ONLY : xpi, xrd, xrv, xg
USE modd_surf_par,   ONLY : xundef
USE modd_forc_atm,  ONLY: xta         ,&! air temperature forcing               (K)
                            xqa       ,&! air specific humidity forcing         (kg/m3)
                            xrhoa     ,&! air density forcing                   (kg/m3)
                            xzs       ,&! orography                             (m)
                            xu        ,&! zonal wind                            (m/s)
                            xv        ,&! meridian wind                         (m/s)
                            xdir_sw   ,&! direct  solar radiation (on horizontal surf.)
                            xsca_sw   ,&! diffuse solar radiation (on horizontal surf.)
                            xlw       ,&! longwave radiation (on horizontal surf.)
                            xps       ,&! pressure at atmospheric model surface (Pa)
                            xpa       ,&! pressure at forcing level             (Pa)
                            xrhoa     ,&! density at forcing level              (kg/m3)
                            xco2      ,&! CO2 concentration in the air          (kg/kg)
                            xsnow     ,&! snow precipitation                    (kg/m2/s)
                            xrain     ,&! liquid precipitation                  (kg/m2/s)
                            xzref       ! height of T,q forcing                 (m)  
USE modd_surfex_omp, ONLY : nindx1sfx,nindx2sfx,nblock,nblocktot, init_dim, reset_dim
!
USE modi_get_luout
USE modi_abor1_sfx
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
#ifdef AIX64
!$ USE OMP_LIB
#endif
!
IMPLICIT NONE
!
!
#ifndef AIX64
!$ INCLUDE 'omp_lib.h'
#endif
!
! global variables
INTEGER,INTENT(IN) :: ksurf_step, knb_atm
INTEGER, DIMENSION(:), INTENT(IN) :: ksize_omp
REAL, DIMENSION(:),INTENT(IN) :: pta1,pta2,pqa1,pqa2,pwind1,pwind2
REAL, DIMENSION(:),INTENT(IN) :: pdir_sw1,pdir_sw2,psca_sw1,psca_sw2,plw1,plw2
REAL, DIMENSION(:),INTENT(IN) :: psnow2,prain2,pps1,pps2,pco21,pco22,pdir1,pdir2

! local variables
REAL :: zdta, zdqa, zddir_sw, zdsca_sw, zdlw,  &
        zdps, zdco2, zdu, zdv, zu1, zv1, zu2, zv2 
REAL :: zpi, znb_atm, zsurf_step, zcoef
INTEGER :: j, inkproma
INTEGER :: iluout
REAL(KIND=JPRB) :: zhook_handle
!========================================================================
!
IF (lhook) CALL dr_hook('OL_TIME_INTERP_ATM',0,zhook_handle)
!
 CALL get_luout('OFFLIN',iluout)
!
zpi = xpi/180.
znb_atm = knb_atm*1.
zsurf_step = ksurf_step*1.-1.
zcoef = zsurf_step / znb_atm
!
!$OMP PARALLEL PRIVATE(INKPROMA,J,ZU1,ZU2,ZV1,ZV2,ZDU,ZDV,ZDTA, &
!$OMP ZDQA,ZDLW,ZDPS,ZDCO2,ZDDIR_SW,ZDSCA_SW)
!
!$ NBLOCK = OMP_GET_THREAD_NUM()
!
IF (nblock==nblocktot) THEN
  CALL init_dim(ksize_omp,0,inkproma,nindx1sfx,nindx2sfx)
ELSE
  CALL init_dim(ksize_omp,nblock,inkproma,nindx1sfx,nindx2sfx)
ENDIF
!
DO j = nindx1sfx,nindx2sfx
  !
  IF (pta1(j)/=xundef) THEN
    ! 
    ! Compute wind components
    !
    ! zonal wind
    zu1 = pwind1(j) * sin(pdir1(j)*zpi)
    zu2 = pwind2(j) * sin(pdir2(j)*zpi)
    zdu = (zu2-zu1)*zcoef
    xu(j) = zu1 + zdu
    !
    zv1 = pwind1(j) * cos(pdir1(j)*zpi) 
    zv2 = pwind2(j) * cos(pdir2(j)*zpi)
    zdv = (zv2-zv1)*zcoef
    xv(j) = zv1 + zdv
    !    
    ! Compute variation from atmospheric time step J and J+1
    !
    zdta     = (pta2(j)-pta1(j))*zcoef
    xta(j) = pta1(j) + zdta
    !
    zdqa     = (pqa2(j)-pqa1(j))*zcoef
    xqa(j) = pqa1(j) + zdqa
    !
    zdlw     = (plw2(j)-plw1(j))*zcoef
    xlw(j) = plw1(j) + zdlw
    !
    zdps     = (pps2(j)-pps1(j))*zcoef
    xps(j) = pps1(j) + zdps
    !
    zdco2    = (pco22(j)-pco21(j))*zcoef
    xco2(j) = pco21(j) + zdco2
    !
    zddir_sw = (pdir_sw2(j)-pdir_sw1(j))*zcoef
    xdir_sw(j,1) = pdir_sw1(j)+zddir_sw
    !
    zdsca_sw = (psca_sw2(j)-psca_sw1(j))*zcoef
    xsca_sw(j,1) = psca_sw1(j)+zdsca_sw
    !
    !
    xrain(j)= prain2(j)
    xsnow(j)= psnow2(j)
    !
    !
    xrhoa(j) = xps(j) / ( xta(j)*xrd * ( 1.+((xrv/xrd)-1.)*xqa(j) ) + xzref(j)*xg )
    !
    ! humidity in kg/m3
    xqa(j) = xqa(j) * xrhoa(j)
    !
  ENDIF
  !
ENDDO
!
 CALL reset_dim(SIZE(pta1),inkproma,nindx1sfx,nindx2sfx)
!
!$OMP END PARALLEL
!
!
! air density
!
! Check No value data
!---------------------
! Error cases
!
IF ((minval(xta)  .EQ.xundef).OR.(minval(xqa).EQ.xundef).OR.&
      (minval(xu).EQ.xundef).OR.(minval(xrain).EQ.xundef).OR.&
      (minval(xsnow).EQ.xundef)) THEN  
    WRITE(iluout,*)'MINVAL(XTA),MINVAL(XQA),MINVAL(XU),MINVAL(XRAIN),MINVAL(XSNOW)'
    WRITE(iluout,*)minval(xta),minval(xqa),minval(xu),minval(xrain),minval(xsnow)
    CALL abor1_sfx('OL_TIME_INTERP_ATM: UNDEFINED VALUE IN ATMOSPHERIC FORCING')
ENDIF
!
IF ((minval(xdir_sw).EQ.xundef).AND.(minval(xsca_sw).EQ.xundef)) THEN
    WRITE(iluout,*)'MINVAL(XSCA_SW),MINVAL(XDIR_SW)'
    WRITE(iluout,*)minval(xsca_sw),minval(xdir_sw)
    CALL abor1_sfx('OL_TIME_INTERP_ATM: UNDEFINED VALUE IN ATMOSPHERIC FORCING')
ENDIF
!
IF ((minval(xps).EQ.xundef).AND.(minval(xzs).EQ.xundef)) THEN
    WRITE(iluout,*)'MINVAL(XPS),MINVAL(XZS)'
    WRITE(iluout,*)minval(xps),minval(xzs)
    CALL abor1_sfx('OL_TIME_INTERP_ATM: UNDEFINED VALUE IN ATMOSPHERIC FORCING')
ENDIF
!
IF (minval(xdir_sw).EQ.xundef) xdir_sw(:,:)=0. ! No direct solar radiation
IF (minval(xsca_sw).EQ.xundef) xsca_sw(:,:)=0. ! No diffuse solar radiation
IF (minval(xps)    .EQ.xundef) THEN            ! No surface Pressure 
   WRITE(iluout,*)' OL_TIME_INTERP_ATM: SURFACE PRESSURE COMPUTED FROM ZS'
   xps(:)  = 101325*(1-0.0065 * xzs(:)/288.15)**5.31
ENDIF
!
!* forcing level pressure from hydrostatism
WHERE(xps(:)/=xundef)
  xpa(:) = xps(:) - xrhoa(:) * xzref(:) * xg
ENDWHERE
!
IF (lhook) CALL dr_hook('OL_TIME_INTERP_ATM',1,zhook_handle)
!
END SUBROUTINE ol_time_interp_atm