make_mask_topd_to_isba.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 make_mask_topd_to_isba (UG, &
                                       ki)
!     #######################
!
!!****  *MAKE_MASK_TOPD_TO_ISBA(*  
!!
!!    PURPOSE
!!    -------
!
!     Create a mask for each catchment. 
!         
!     
!!**  METHOD
!!    ------
!
!!    EXTERNAL
!!    --------
!!
!!    none
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------ 
!!                     
!!    REFERENCE
!!    ---------
     
!!    AUTHOR
!!    ------
!!
!!      L. Labatut & K. Chancibault     * CNRM *
!!
!!    MODIFICATIONS
!!    -------------
!!
!!      Original   16/03/2005
!!                 03/2014 (E. Artinian) manages the option CGRID='IGN'
!!                 07/2015 (E. Artinian) corrections  for option CGRID='IGN'
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
!
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
!
USE modd_topodyn,       ONLY: nncat, nnyc, xy0, xdxt, nnxc,&
                                xx0, xtopd, xnul, nline, nmesht
USE modd_coupling_topd, ONLY: nmaskt
USE modd_surf_par,        ONLY : xundef, nundef
!
USE modi_get_luout
USE modi_abor1_sfx
USE modi_write_file_map
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
TYPE(surf_atm_grid_t), INTENT(INOUT) :: ug
!
INTEGER, INTENT(IN) :: ki    ! Grid dimensions
!
!*      0.2    declarations of local variables
!
 CHARACTER(LEN=30)  :: yvar        ! name of results file
INTEGER            :: jcat, jj, ji, idx     ! loop control  
INTEGER            :: ii,iline ! work integer variables
INTEGER            :: idxm     ! indexes of Isba grid meshes and nodes
INTEGER            :: iluout   ! unit
REAL               :: zxt, zyt ! catchment grid nodes Lambert II coordinates 
REAL               :: zx1, zx2, zx3, zx4, zy1, zy2, zy3, zy4  ! Isba mesh Lambert II coordinates
REAL               :: zxa, zxb, zya, zyb
REAL, DIMENSION(NNCAT,NMESHT):: zwrk
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('MAKE_MASK_TOPD_TO_ISBA',0,zhook_handle)
!
 CALL get_luout('OFFLIN',iluout)
!  
!WRITE(*,*) 'Y0=', XY0(1), XY0(2)
!WRITE(*,*) 'X0=', XX0(1), XX0(2)
!WRITE(*,*) 'X et Y isba=', XXI(1), XXI(180), XYI(1), XYI(180)
!loop on catchments
DO jcat=1,nncat
  !
  idxm  = 1
  CALL init_4points(idxm,zx1,zx2,zx3,zx4,zy1,zy2,zy3,zy4)
  !
  DO jj=1,nnyc(jcat)   ! number of topographic points on the y axis
    !
    zyt = xy0(jcat) + (jj-1) * xdxt(jcat)
    zyt = zyt + 0.5 * xdxt(jcat)
    !
    DO ji=1,nnxc(jcat)
      !
      zxt = xx0(jcat) + (ji-1) * xdxt(jcat)
      zxt = zxt + 0.5 * xdxt(jcat)
      !
      idx = (jj-1) * nnxc(jcat) + ji   ! index of the point among all in the catchment
      !
      !* on vérifie que le pixel du MNT appartient au BV
      IF ( xtopd(jcat,idx).NE.xnul(jcat) ) THEN
        !* calcule des coordonées X et Y inf et sup de la maille ISBA considérée
        CALL get_coord(zxt,zyt,zx1,zx2,zx3,zx4,zy1,zy2,zy3,zy4,zxa,zya,zxb,zyb)
        !* si on se trouve sur le premier pixel du MNT ou si le pixel du MNT n'est pas 
        !dans la maille Isba considérée (qui est celle dans laquelle se trouve le pixel précédent)
        IF (zxt.LT.zxa.OR.zxt.GE.zxb.OR.zyt.LT.zya.OR.zyt.GE.zyb) THEN
          !* on repart de la première maille de la grille Isba
          idxm = 1
          CALL init_4points(idxm,zx1,zx2,zx3,zx4,zy1,zy2,zy3,zy4)
          CALL get_coord(zxt,zyt,zx1,zx2,zx3,zx4,zy1,zy2,zy3,zy4,zxa,zya,zxb,zyb)
          !* on parcours les mailles de la grille Isba, jusqu'à ce qu'on trouve la maille à laquelle appartient le pixel du MNT
          DO WHILE (zxt.LT.zxa.OR.zxt.GE.zxb.OR.zyt.LT.zya.OR.zyt.GE.zyb)
            idxm = idxm + 1
            IF (idxm.GE.ki) THEN
              WRITE(*,*) 'ZXT', zxt,'ZYT',zyt
              WRITE(*,*) 'indices Isba:',idxm,'>=',ki
              CALL abor1_sfx("MAKE_MASK_TOPD_TO_ISBA: PROBLEM")
            ENDIF
            CALL init_4points(idxm,zx1,zx2,zx3,zx4,zy1,zy2,zy3,zy4)
            CALL get_coord(zxt,zyt,zx1,zx2,zx3,zx4,zy1,zy2,zy3,zy4,zxa,zya,zxb,zyb)
          ENDDO
        ENDIF
        IF (nline(jcat,idx)/=0) nmaskt(jcat,nline(jcat,idx)) = idxm
      ENDIF
    ENDDO
  ENDDO
ENDDO
!
yvar='.mask_topd'
WHERE (nmaskt(:,:)/=nundef)
  zwrk(:,:)=REAL(nmaskt)
ELSEWHERE
  zwrk(:,:)=xundef
ENDWHERE
 CALL write_file_map(zwrk,yvar)
!CALL WRITE_FILE_MAP(REAL(NMASKT),YVAR)
!
IF (lhook) CALL dr_hook('MAKE_MASK_TOPD_TO_ISBA',1,zhook_handle)
!
 CONTAINS
!
SUBROUTINE init_4points(KDXM,PX1,PX2,PX3,PX4,PY1,PY2,PY3,PY4)
!
USE modd_coupling_topd, ONLY: nimax, xxi, xyi
USE mode_gridtype_ign
!
INTEGER, INTENT(IN) :: kdxm
REAL, INTENT(OUT) :: px1, px2, px3, px4
REAL, INTENT(OUT) :: py1, py2, py3, py4
REAL, DIMENSION(KI)    :: zdx, zdy
!
INTEGER :: iline, ii, idxn
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('MAKE_MASK_TOPD_TO_ISBA:INIT_4POINTS',0,zhook_handle)
!
IF (ug%CGRID=='IGN') THEN 
  CALL get_gridtype_ign(ug%XGRID_PAR,pdx=zdx,pdy=zdy)
  idxn=kdxm
  !on va juste retourner les quatre coins de la maille, les XXI et XYI etant les coordonees du centre
  !on se contente de verifier si la maille TOP est dans la maille SURFEX
  !la grille est reguliere dans les deux sens mais pas forcement rectangulaire
  px1=xxi(idxn)-zdx(idxn)/2.0
  px2=xxi(idxn)+zdx(idxn)/2.0
  px3=xxi(idxn)-zdx(idxn)/2.0
  px4=xxi(idxn)+zdx(idxn)/2.0
  py1=xyi(idxn)-zdy(idxn)/2.0
  py2=xyi(idxn)-zdy(idxn)/2.0
  py3=xyi(idxn)+zdy(idxn)/2.0
  py4=xyi(idxn)+zdy(idxn)/2.0
ELSE
 iline = int(kdxm/(nimax))+1      ! number of the current line
 ii    = kdxm-((iline-1)*nimax)   ! index of point in the line
 idxn  = (iline-1)*(nimax+1)+ii   ! indice du point dans la grille isba
!
 px1 = xxi(idxn)              ! coordonnée X du point courant
 px2 = xxi(idxn+1)            ! coordonnée X du point suivant
 px3 = xxi(idxn+(nimax+1))    ! coordonnée X du point aligné sur la ligne suivante 
 px4 = xxi(idxn+1+(nimax+1))  ! coordonnée X du point suivant sur la ligne suivante
!
 py1 = xyi(idxn)              ! coordonnée Y du point courant
 py2 = xyi(idxn+1)            ! coordonnée Y du point suivant
 py3 = xyi(idxn+(nimax+1))    ! coordonnée Y du point aligné sur la ligne suivante
 py4 = xyi(idxn+1+(nimax+1))  ! coordonnée Y du point suivant sur la ligne suivante
ENDIF
!
IF (lhook) CALL dr_hook('MAKE_MASK_TOPD_TO_ISBA:INIT_4POINTS',1,zhook_handle)
!
END SUBROUTINE init_4points
!
SUBROUTINE get_coord(PXT,PYT,PX1,PX2,PX3,PX4,PY1,PY2,PY3,PY4,&
                     pxa,pya,pxb,pyb)
!
REAL, INTENT(IN) :: pxt, pyt
REAL, INTENT(IN) :: px1, px2, px3, px4
REAL, INTENT(IN) :: py1, py2, py3, py4
REAL, INTENT(OUT) :: pxa, pya, pxb, pyb
REAL :: zfa, zfb
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('MAKE_MASK_TOPD_TO_ISBA:GET_COORD',0,zhook_handle)
!
IF((px3-px1).EQ.0.0) THEN
  pxa = px1
ELSE
  CALL get_line_param(px1,py1,px3,py3,zfa,zfb)
  pxa = (zyt-zfb)/zfa
ENDIF
!
IF ((px4-px2).EQ.0.0) THEN
  pxb = px2
ELSE
  CALL get_line_param(px2,py2,px4,py4,zfa,zfb)
  pxb = (zyt-zfb)/zfa
ENDIF
!
IF ((py2-py1).EQ.0.0) THEN
  pya = py2
ELSE
  CALL get_line_param(px1,py1,px2,py2,zfa,zfb)
  pya = zfa*zxt+zfb
ENDIF
!
IF ((py4-py3).EQ.0.0) THEN
  pyb = py4
ELSE
  CALL get_line_param(px3,py3,px4,py4,zfa,zfb)
  pyb = zfa*zxt+zfb
ENDIF
!
IF (lhook) CALL dr_hook('MAKE_MASK_TOPD_TO_ISBA:GET_COORD',1,zhook_handle)
!
END SUBROUTINE get_coord
!
SUBROUTINE get_line_param(PX1,PY1,PX2,PY2,PFA,PFB)
!
REAL, INTENT(IN) :: px1, px2, py1, py2
REAL, INTENT (OUT) :: pfa, pfb
REAL(KIND=JPRB) :: zhook_handle
!
IF (lhook) CALL dr_hook('MAKE_MASK_TOPD_TO_ISBA:GET_LINE_PARAM',0,zhook_handle)
!
pfa = (py2 - py1) / (px2 - px1)  ! slope
pfb = py1 - pfa * px1 ! offset
!
IF (lhook) CALL dr_hook('MAKE_MASK_TOPD_TO_ISBA:GET_LINE_PARAM',1,zhook_handle)
!
END SUBROUTINE get_line_param
!
END SUBROUTINE make_mask_topd_to_isba