subscale_aos.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.
!     #########
      SUBROUTINE subscale_aos (UG, USS, &
                               oz0effi,oz0effj,psea)
!     #############################################
!
!!*SUBSCALE_AOS  computes the sum of the ratio: (h'-h)/L when  h'/L >h/L  
!!                  the ' is for subgrid scale orography
!!
!!
!!    METHOD
!!    ------
!!    See M.Georgelin and al. July 1994, Monthly Weather Review.
!!   
!!    EXTERNAL
!!    --------
!!
!!    AUTHOR
!!    ------
!!
!!    M. Georgelin      Laboratoire d'Aerologie
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    18/12/95
!!    (V. Masson) 10/03/97 rewrites the routine in doctor norm.
!!                         computations are made even if a only a few subsquares
!!                         contains data points.
!!                         returns to the calling routine the localization of
!!                         the points where the z0 coefficients are available.
!!    (V. Masson) 03/2004  externalization
!!
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
!
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_surf_atm_sso_n, ONLY : surf_atm_sso_t
!
USE modd_pgdwork,        ONLY : nsso, xssqo, lssqo
USE modd_pgd_grid,       ONLY : nl, cgrid, xgrid_par, ngrid_par
!
USE modi_get_adjacent_meshes
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_get_mesh_dim
IMPLICIT NONE
!
!*    0.1    Declaration of dummy arguments
!            ------------------------------
!
!
TYPE(surf_atm_grid_t), INTENT(INOUT) :: ug
TYPE(surf_atm_sso_t), INTENT(INOUT) :: uss
!
LOGICAL, DIMENSION(:), INTENT(OUT) :: oz0effi! .T. : the z0eff coefficients
!                                            ! are computed at grid point
!                                            ! .F. : not enough sub-grid
!                                            ! information avalaible to
!                                            ! compute the coefficients
LOGICAL, DIMENSION(:), INTENT(OUT) :: oz0effj! .T. : the z0eff coefficients
!                                            ! are computed at grid point
!                                            ! .F. : not enough sub-grid
!                                            ! information avalaible to
!                                            ! compute the coefficients
REAL,    DIMENSION(:), INTENT(IN)  :: psea   ! sea fraction

!
!*    0.2    Declaration of indexes
!            ----------------------
!
!
INTEGER :: jl          ! loop index on grid meshs
INTEGER :: il          ! grid mesh index of second subgrid point used
INTEGER :: jiss, jjss  ! loop indexes for subsquares arrays
INTEGER :: jnext       ! loop index on subgrid meshes
INTEGER :: inext       ! index to add to JISS or JJSS to obtain the following
!                      ! point containing a data in a segment
INTEGER, DIMENSION(NL) :: ileft   ! index of left   grid mesh 
INTEGER, DIMENSION(NL) :: iright  ! index of right  grid mesh 
INTEGER, DIMENSION(NL) :: itop    ! index of top    grid mesh 
INTEGER, DIMENSION(NL) :: ibottom ! index of bottom grid mesh
!
!*    0.3    Declaration of counters inside a grid (JL)
!            -----------------------
!
INTEGER :: iho2counterjp ! number of times where h/2 has been summed for JP coef
INTEGER :: iho2counterjm ! number of times where h/2 has been summed for JM coef
INTEGER :: iho2counterip ! number of times where h/2 has been summed for IP coef
INTEGER :: iho2counterim ! number of times where h/2 has been summed for IM coef
INTEGER :: iaoscounter   ! number of segments where A/S has been summed
INTEGER :: iaosdist ! distance between first and last subsquares used in
!                   ! computation of A/S in a subsegment of the grid
LOGICAL :: gfirst   ! T indicates the first point has been found for this segment.
!
!*    0.4    Declaration of working arrays inside a grid (JL)
!            -----------------------------
!
REAL, DIMENSION(NSSO) :: zaosip ! A/S in each subsegment for IP coef.
REAL, DIMENSION(NSSO) :: zaosim ! A/S in each subsegment for IM coef.
REAL, DIMENSION(NSSO) :: zaosjp ! A/S in each subsegment for JP coef.
REAL, DIMENSION(NSSO) :: zaosjm ! A/S in each subsegment for JM coef.
REAL :: zaip      ! Area in a subsegment for IP coef.
REAL :: zaim      ! Area in a subsegment for IM coef.
REAL :: zajp      ! Area in a subsegment for JP coef.
REAL :: zajm      ! Area in a subsegment for JM coef.
REAL :: zsumho2ip ! sum of h/2 in the grid for IP coef.
REAL :: zsumho2im ! sum of h/2 in the grid for IM coef.
REAL :: zsumho2jp ! sum of h/2 in the grid for JP coef.
REAL :: zsumho2jm ! sum of h/2 in the grid for JM coef.
REAL :: zssaos    ! A/S between 2 following points along a subsegment
REAL :: zslope    ! slope between 2 following points along a subsegment
REAL :: zdxeff    ! width of a subsquare along I axis
REAL :: zdyeff    ! width of a subsquare along J axis
!
!*    0.5    Declaration of other local variables
!            ------------------------------------
!
REAL, DIMENSION(NL)   :: zdx      ! grid mesh size in x direction
REAL, DIMENSION(NL)   :: zdy      ! grid mesh size in y direction
REAL, DIMENSION(0:NL) :: zslopeip ! x mean slope
REAL, DIMENSION(0:NL) :: zslopejp ! y mean slope
REAL(KIND=JPRB) :: zhook_handle
!----------------------------------------------------------------------------
!
!*    1.     Initializations
!            ---------------
!
!*    1.1    Occurence of computation of the coefficients
!            --------------------------------------------
!
IF (lhook) CALL dr_hook('SUBSCALE_AOS',0,zhook_handle)
oz0effi(:)=.false.
oz0effj(:)=.false.
!
!*    1.2    Grid dimension (meters)
!            -----------------------
!
 CALL get_mesh_dim(cgrid,ngrid_par,nl,xgrid_par,zdx,zdy,ug%XMESH_SIZE)
!
!
!*    1.3    Left, top, right and bottom adjacent gris meshes
!            ------------------------------------------------
!
 CALL get_adjacent_meshes(cgrid,ngrid_par,nl,xgrid_par,ileft,iright,itop,ibottom)
!
!
!*    1.4    Mean slopes between 2 grid meshes
!            -----------
!
zslopeip(0) = 0.
zslopejp(0) = 0.
!
DO jl=1,nl
  IF (iright(jl)/=0 .AND. ileft(jl)/=0) THEN
    zslopeip(jl) =  0.5 * ( uss%XAVG_ZS(iright(jl)) - uss%XAVG_ZS(jl) ) &
                          / ( 0.5 * (zdx(iright(jl)) + zdx(jl)) ) &
                    + 0.5 * ( uss%XAVG_ZS(jl) - uss%XAVG_ZS(ileft(jl)) ) &
                          / ( 0.5 * (zdx(jl)  + zdx(ileft(jl))) )  
  ELSE
    zslopeip(jl) = 0.
  END IF
  IF (itop(jl)/=0 .AND. ibottom(jl)/=0) THEN
    zslopejp(jl) =  0.5 * ( uss%XAVG_ZS(itop(jl))     - uss%XAVG_ZS(jl) ) &
                          / ( 0.5 * (zdy(itop(jl))     + zdy(jl)) ) &
                    + 0.5 * ( uss%XAVG_ZS(jl) - uss%XAVG_ZS(ibottom(jl)) ) &
                          / ( 0.5 * (zdy(jl)  + zdy(ibottom(jl))) )  
  ELSE
    zslopejp(jl) = 0.
  END IF
END DO
!
!----------------------------------------------------------------------------
!
!*    2.     Loop on grid points
!            -------------------
!
DO jl=1,nl
!
!*    2.1    No land in grid mesh
!            --------------------
!
    IF (psea(jl)==1.) cycle
!
!*    2.2    Index Initializations
!            ---------------------
!
    zdxeff=zdx(jl)/float(nsso)
    zdyeff=zdy(jl)/float(nsso)
!
!----------------------------------------------------------------------------
!
!*    3.     Computations for IP and IM fields
!            ---------------------------------
!
    zaosip(:)=0.
    zaosim(:)=0.
    zsumho2ip=0.
    zsumho2im=0.
    iho2counterip=0
    iho2counterim=0
    iaoscounter=0
!
!*    3.1    loop on jss index (there is no specific computation along j)
!            -----------------
!
    DO jjss=1,nsso
!
!*    3.1.1  initializes counters for the A/S subscale segment computation
!
      gfirst = .true.
      iaosdist=0
      zaip=0.
      zaim=0.
!
!*    3.2    loop on iss index
!            -----------------
!
      DO jiss=1,nsso
!
!*    3.3    search for two consecutive grid points
!            --------------------------------------
!
!*    3.3.1 first one
!
        IF (.NOT. lssqo(jiss,jjss,jl) ) cycle
!
!*    3.3.2  second one (up to one grid mesh further)
!
        DO jnext=1,nsso
          IF (jiss+jnext>nsso) THEN
            il    = iright(jl)
            inext = jiss+jnext-nsso
          ELSE
            il    = jl
            inext = jiss+jnext
          END IF
          ! no right point
          IF (il==0) EXIT
          ! subgrid data found
          IF (lssqo(inext,jjss,il)) EXIT
        END DO
!
!*    3.3.3  none found: end of loop along jss
!                        ---------------------
!
        IF (jnext>=nsso+1) EXIT
!
!*    3.3.4  second point outside of the domain: end of loop along iss
!                                                ---------------------
!
        IF (il==0) EXIT
!
!*    3.4    add terms to sums of A/S and h/2
!            --------------------------------
!
        IF (gfirst) iaoscounter=iaoscounter+1
        gfirst = .false.
        iaosdist   =iaosdist+jnext
!
!*    3.4.1  mean slope
!
        zslope=zslopeip(jl)
!
!*    3.4.2  A/S term
!
        zssaos =  xssqo(inext,jjss,il) - xssqo(jiss,jjss,jl) &
                  - zslope * zdxeff * jnext  
        IF (zssaos>0.) zaip=zaip+zssaos
        IF (zssaos<0.) zaim=zaim-zssaos
!
!*    3.4.3  h/2 term
!
        IF (zssaos>0.) THEN
          zsumho2ip = zsumho2ip + 0.5 * zssaos
          iho2counterip=iho2counterip+1
        END IF
        IF (zssaos<0.) THEN
          zsumho2im = zsumho2im - 0.5 * zssaos
          iho2counterim=iho2counterim+1
        END IF
!
!*    3.5    end of loop on iss index
!            ------------------------
!
      END DO
      IF (iaosdist>0) THEN
        zaosip(jjss)=zaip/(zdxeff*iaosdist)
        zaosim(jjss)=zaim/(zdxeff*iaosdist)
      END IF
!
!*    3.6    end of loop on jss index
!            ------------------------
!
    END DO
!
!*    3.7    end of IP and IM coefficients
!            -----------------------------
!
    IF (iaoscounter>0) THEN
      uss%XAOSIP(jl)=sum(zaosip) / iaoscounter
      uss%XAOSIM(jl)=sum(zaosim) / iaoscounter
      IF (iho2counterip>0) THEN
        uss%XHO2IP(jl)=zsumho2ip   / iho2counterip
      ELSE
        uss%XHO2IP(jl)=0.
      END IF
      IF (iho2counterim>0) THEN
        uss%XHO2IM(jl)=zsumho2im   / iho2counterim
      ELSE
        uss%XHO2IM(jl)=0.
      END IF
      oz0effi(jl)=.true.
    END IF
!
!----------------------------------------------------------------------------
!
!*    4.     Computations for JP and JM fields
!            ---------------------------------
!
    zaosjp(:)=0.
    zaosjm(:)=0.
    zsumho2jp=0.
    zsumho2jm=0.
    iho2counterjp=0
    iho2counterjm=0
    iaoscounter=0
!
!*    4.1    loop on iss index (there is no specific computation along i)
!            -----------------
!
    DO jiss=1,nsso
!
!*    4.1.1  initializes counters for the A/S subscale segment computation
!
      gfirst = .true.
      iaosdist=0
      zajp=0.
      zajm=0.
!
!*    4.2    loop on jss index
!            -----------------
!
      DO jjss=1,nsso
!
!*    4.3    search for two consecutive grid points
!            --------------------------------------
!
!*    4.3.1 first one
!
        IF (.NOT. lssqo(jiss,jjss,jl) ) cycle
!
!*    4.3.2  second one (up to one grid mesh further)
!
        DO jnext=1,nsso
          IF (jjss+jnext>nsso) THEN
            il    = itop(jl)
            inext = jjss+jnext-nsso
          ELSE
            il    = jl
            inext = jjss+jnext
          END IF
          ! no right point
          IF (il==0) EXIT
          ! subgrid data found
          IF (lssqo(jiss,inext,il)) EXIT
        END DO
!
!*    4.3.3  none found: end of loop along jss
!                        ---------------------
!
        IF (jnext>=nsso+1) EXIT
!
!*    4.3.4  second point outside of the domain: end of loop along iss
!                                                ---------------------
!
        IF (il==0) EXIT
!
!
!*    4.4    add terms to sums of A/S and h/2
!            --------------------------------
!
        IF (gfirst) iaoscounter=iaoscounter+1
        gfirst = .false.
        iaosdist   =iaosdist+jnext
!
!*    4.4.1  mean slope
!
        zslope=zslopejp(jl)
!
!*    4.4.2  A/S term
!
        zssaos =  xssqo(jiss,inext,il) - xssqo(jiss,jjss,jl) &
                  - zslope * zdyeff * jnext  
        IF (zssaos>0.) zajp=zajp+zssaos
        IF (zssaos<0.) zajm=zajm-zssaos
!
!*    4.4.3  h/2 term
!
        IF (zssaos>0.) THEN
          zsumho2jp = zsumho2jp + 0.5 * zssaos
          iho2counterjp=iho2counterjp+1
        END IF
        IF (zssaos<0.) THEN
          zsumho2jm = zsumho2jm - 0.5 * zssaos
          iho2counterjm=iho2counterjm+1
        END IF
!
!*    4.5    end of loop on jss index
!            ------------------------
!
      END DO
      IF (iaosdist>0) THEN
        zaosjp(jiss)=zajp/(zdyeff*iaosdist)
        zaosjm(jiss)=zajm/(zdyeff*iaosdist)
      END IF
!
!*    4.6    end of loop on iss index
!            ------------------------
!
    END DO
!
!*    4.7    end of JP and JM coefficients
!            -----------------------------
!
    IF (iaoscounter>0) THEN
      uss%XAOSJP(jl)=sum(zaosjp) /iaoscounter
      uss%XAOSJM(jl)=sum(zaosjm) /iaoscounter
      IF (iho2counterjp>0) THEN
        uss%XHO2JP(jl)=zsumho2jp   /iho2counterjp
      ELSE
        uss%XHO2JP(jl)=0.
      END IF
      IF (iho2counterjm>0) THEN
        uss%XHO2JM(jl)=zsumho2jm   /iho2counterjm
      ELSE
        uss%XHO2JM(jl)=0.
      END IF
      oz0effj(jl)=.true.
    END IF
!
!-------------------------------------------------------------------------------
!
!*    5.     Next grid point
!            ---------------
!
END DO
IF (lhook) CALL dr_hook('SUBSCALE_AOS',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE subscale_aos