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 (U, UG, USS, OZ0EFFI, OZ0EFFJ)
!     #############################################
!
!!*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_n, ONLY : surf_atm_t
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
USE modd_sso_n, ONLY : sso_t
!
USE modd_surfex_mpi, ONLY : nrank, npio
USE modd_pgdwork,        ONLY : nsso, xssqo, lssqo
USE modd_pgd_grid,       ONLY : nl
!
USE modi_get_adjacent_meshes
USE modi_read_and_send_mpi
USE modi_gather_and_write_mpi
!
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_t), INTENT(INOUT) :: U
TYPE(surf_atm_grid_t), INTENT(INOUT) :: UG
TYPE(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
!
!*    0.2    Declaration of indexes
!            ----------------------
!
INTEGER :: JS1, JS2
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(:), ALLOCATABLE :: ILEFT   ! index of left   grid mesh 
INTEGER, DIMENSION(:), ALLOCATABLE :: IRIGHT  ! index of right  grid mesh 
INTEGER, DIMENSION(:), ALLOCATABLE :: ITOP    ! index of top    grid mesh 
INTEGER, DIMENSION(:), ALLOCATABLE :: 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
!            ------------------------------------
!
INTEGER, DIMENSION(:), ALLOCATABLE :: ISSO, ISSOT
INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: ISSQOT, ISSQO
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZSSQO
LOGICAL, DIMENSION(:), ALLOCATABLE :: GZ0EFFI, GZ0EFFJ
REAL, DIMENSION(:), ALLOCATABLE   :: ZDX, ZDY, ZSEA, ZAVG_ZS, ZMESH_SIZE    ! grid mesh size in x direction
REAL, DIMENSION(:), ALLOCATABLE   :: ZAOSIP_ALL, ZAOSIM_ALL, ZAOSJP_ALL, ZAOSJM_ALL 
REAL, DIMENSION(:), ALLOCATABLE   :: ZHO2IP_ALL, ZHO2IM_ALL, ZHO2JP_ALL, ZHO2JM_ALL 
REAL :: ZSLOPEIP ! x mean slope
REAL :: 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.
!
IF (nrank==npio) THEN
  ALLOCATE(zsea(u%NDIM_FULL))
  ALLOCATE(zavg_zs(u%NDIM_FULL))
  ALLOCATE(zmesh_size(u%NDIM_FULL))  
  ALLOCATE(issqot(u%NDIM_FULL,nsso,nsso))
  ALLOCATE(zssqo(u%NDIM_FULL,nsso,nsso))
ELSE
  ALLOCATE(zsea(0))
  ALLOCATE(zavg_zs(0)) 
  ALLOCATE(issqot(0,0,0))
  ALLOCATE(zssqo(0,0,0))
ENDIF
!
 CALL gather_and_write_mpi(u%XSEA,zsea)
 CALL gather_and_write_mpi(uss%XAVG_ZS,zavg_zs)
 CALL gather_and_write_mpi(ug%G%XMESH_SIZE,zmesh_size)  
 CALL gather_and_write_mpi(xssqo,zssqo)
!
ALLOCATE(issqo(SIZE(xssqo,1),nsso,nsso))
issqo(:,:,:) = 0
WHERE (lssqo(:,:,:)) issqo(:,:,:) = 1
 CALL gather_and_write_mpi(issqo,issqot)
DEALLOCATE(issqo)
!
!*    1.2    Grid dimension (meters)
!            -----------------------
!
IF (nrank==npio) THEN
  !
  ALLOCATE(zdx(u%NDIM_FULL),zdy(u%NDIM_FULL))
  CALL get_mesh_dim(ug%G%CGRID,ug%NGRID_FULL_PAR,u%NDIM_FULL,ug%XGRID_FULL_PAR,zdx,zdy,zmesh_size)
  DEALLOCATE(zmesh_size)
  !
  !
  !*    1.3    Left, top, right and bottom adjacent gris meshes
  !            ------------------------------------------------
  !
  ALLOCATE(ileft(u%NDIM_FULL),iright(u%NDIM_FULL),itop(u%NDIM_FULL),ibottom(u%NDIM_FULL))
  CALL get_adjacent_meshes(ug%G%CGRID,ug%NGRID_FULL_PAR,u%NDIM_FULL,ug%XGRID_FULL_PAR,&
                                ileft,iright,itop,ibottom)
  !
  ALLOCATE(gz0effi(u%NDIM_FULL),gz0effj(u%NDIM_FULL))
  ALLOCATE(zho2im_all(u%NDIM_FULL),zho2ip_all(u%NDIM_FULL),&
           zaosim_all(u%NDIM_FULL),zaosip_all(u%NDIM_FULL))
  ALLOCATE(zho2jm_all(u%NDIM_FULL),zho2jp_all(u%NDIM_FULL),&
           zaosjm_all(u%NDIM_FULL),zaosjp_all(u%NDIM_FULL))   
  gz0effi(:) = .false.
  gz0effj(:) = .false.
  !
  !*    1.4    Mean slopes between 2 grid meshes
  !            -----------
  !
  !
  DO jl=1,u%NDIM_FULL
    !
    zslopeip = 0.
    zslopejp = 0.
    !
    IF (iright(jl)/=0 .AND. ileft(jl)/=0) THEN
      zslopeip =  0.5 * ( zavg_zs(iright(jl)) - zavg_zs(jl) ) &
                            / ( 0.5 * (zdx(iright(jl)) + zdx(jl)) ) &
                      + 0.5 * ( zavg_zs(jl) - zavg_zs(ileft(jl)) ) &
                            / ( 0.5 * (zdx(jl)  + zdx(ileft(jl))) )  
    ELSE
      zslopeip = 0.
    END IF
    IF (itop(jl)/=0 .AND. ibottom(jl)/=0) THEN
      zslopejp =  0.5 * ( zavg_zs(itop(jl))     - zavg_zs(jl) ) &
                            / ( 0.5 * (zdy(itop(jl))     + zdy(jl)) ) &
                      + 0.5 * ( zavg_zs(jl) - zavg_zs(ibottom(jl)) ) &
                            / ( 0.5 * (zdy(jl)  + zdy(ibottom(jl))) )  
    ELSE
      zslopejp = 0.
    END IF
!
!----------------------------------------------------------------------------
!
!*    2.     Loop on grid points
!            -------------------
!!
!*    2.1    No land in grid mesh
!            --------------------
!
    IF (zsea(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 (issqot(jl,jiss,jjss)==0 ) 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 (issqot(il,inext,jjss)==1) 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
!
!*    3.4.2  A/S term
!
        zssaos =  zssqo(il,inext,jjss) - zssqo(jl,jiss,jjss) &
                  - 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
      zaosip_all(jl)=sum(zaosip) / iaoscounter
      zaosim_all(jl)=sum(zaosim) / iaoscounter
      IF (iho2counterip>0) THEN
        zho2ip_all(jl)=zsumho2ip   / iho2counterip
      ELSE
        zho2ip_all(jl)=0.
      END IF
      IF (iho2counterim>0) THEN
        zho2im_all(jl)=zsumho2im   / iho2counterim
      ELSE
        zho2im_all(jl)=0.
      END IF
      gz0effi(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 (issqot(jl,jiss,jjss)==0 ) 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 (issqot(il,jiss,inext)==1) 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
!
!*    4.4.2  A/S term
!
        zssaos =  zssqo(il,jiss,inext) - zssqo(jl,jiss,jjss) &
                  - 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
      zaosjp_all(jl)=sum(zaosjp) /iaoscounter
      zaosjm_all(jl)=sum(zaosjm) /iaoscounter
      IF (iho2counterjp>0) THEN
        zho2jp_all(jl)=zsumho2jp   /iho2counterjp
      ELSE
        zho2jp_all(jl)=0.
      END IF
      IF (iho2counterjm>0) THEN
        zho2jm_all(jl)=zsumho2jm   /iho2counterjm
      ELSE
        zho2jm_all(jl)=0.
      END IF
      gz0effj(jl)=.true.
    END IF
!
  END DO
  ! 
ELSE
  ALLOCATE(zho2im_all(0),zho2ip_all(0),zaosim_all(0),zaosip_all(0))
  ALLOCATE(zho2jm_all(0),zho2jp_all(0),zaosjm_all(0),zaosjp_all(0))
ENDIF
!
DEALLOCATE(zssqo,issqot)
DEALLOCATE(zsea)  
DEALLOCATE(zavg_zs)
!
 CALL read_and_send_mpi(zaosip_all,uss%XAOSIP)
 CALL read_and_send_mpi(zaosim_all,uss%XAOSIM)
 CALL read_and_send_mpi(zho2ip_all,uss%XHO2IP)
 CALL read_and_send_mpi(zho2im_all,uss%XHO2IM)
DEALLOCATE(zho2im_all,zho2ip_all,zaosim_all,zaosip_all)
!
 CALL read_and_send_mpi(zaosjp_all,uss%XAOSJP)
 CALL read_and_send_mpi(zaosjm_all,uss%XAOSJM)
 CALL read_and_send_mpi(zho2jp_all,uss%XHO2JP)
 CALL read_and_send_mpi(zho2jm_all,uss%XHO2JM)
DEALLOCATE(zho2jm_all,zho2jp_all,zaosjm_all,zaosjp_all)
!
IF (nrank==npio) THEN
  ALLOCATE(issot(u%NDIM_FULL))
  issot(:) = 0
  WHERE (gz0effi(:)) issot(:) = 1
  DEALLOCATE(gz0effi)
ELSE
  ALLOCATE(issot(0))
ENDIF
ALLOCATE(isso(nl))
 CALL read_and_send_mpi(issot,isso)
WHERE(isso(:)==1) oz0effi(:) = .true.
!
IF (nrank==npio) THEN
  issot(:) = 0
  WHERE (gz0effj(:)) issot(:) = 1
  DEALLOCATE(gz0effj)
ENDIF
 CALL read_and_send_mpi(issot,isso)
WHERE(isso(:)==1) oz0effj(:) = .true.
DEALLOCATE(isso)
!
DEALLOCATE(issot)
!
!-------------------------------------------------------------------------------
!
!*    5.     Next grid point
!            ---------------
!

IF (lhook) CALL dr_hook('SUBSCALE_AOS',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE subscale_aos