MODULE getlandseamask
  !
  !================================================================================================================================
  !! MODULE   : getlandseamask
  !!
  !> BRIEF      This module contains a few methods in order to build a land surface mask which can serve
  !!            to test the routing and its ability to build a river network.
  !!
  !! DESCRIPTION  : None
  !!
  !! RECENT CHANGE(S): None
  !!
  !! MAIN OUTPUT VARIABLE(S):
  !!
  !! REFERENCES   : None
  !!
  !! FLOWCHART    : This module has 2 possible modes for building the land/sea mask :
  !!                TOPOFILE : This configuration option contains a file name of the orography file
  !!                           which will be used in order to build the land/sea mask.
  !!                IIM_LON, JJM_lon : the number of points the grid should have in longitude and
  !!                                   latitude. This is only used when we build the mask based on orography.
  !!                                   Should these parameters not be provided, then the gird of the orography
  !!                                   file is used.
  !!                CONTFRACFILE :: File name in which we will find variable "Contfrac". This land/sea mask will be 
  !!                                used with the associated lat/lon grid.
  !!                WEST_EAST,  SOUTH_NORTH : The zoom to be performed on the grid so as to retain only this window
  !!                                          of the land-sea mask.
  !! 
  !! $HeadURL: svn://forge.ipsl.fr/orchidee/trunk/ORCHIDEE/src_driver/getlandseamask.f90 $
  !! $Date: 2023-07-04 15:40:20 +0200 (mar. 04 juil. 2023) $
  !! $Revision: 8069 $
  !_================================================================================================================================

  USE ioipsl_para
  USE mod_orchidee_para
  USE control
  USE constantes_soil_var
  USE constantes_var
  USE constantes

  USE netcdf

  IMPLICIT NONE

  PRIVATE
  PUBLIC :: getlandseamask_init, getlandseamask_read

  INTEGER(i_std), SAVE                           :: iim, jjm, nbland
  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: lon, lat
  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: contfrac, orog

CONTAINS
  !================================================================================================================================
  !! SUBROUTINE   : getlandseamask_init
  !!
  !> BRIEF         Builds the land/sea mask and returns its size to the calling program.
  !!
  !! DESCRIPTION  : None
  !!
  !! RECENT CHANGE(S): None
  !!
  !! MAIN OUTPUT VARIABLE(S):
  !!
  !! REFERENCES   : None
  !!
  !! FLOWCHART    : None
  !! \n
  !_================================================================================================================================
  SUBROUTINE getlandseamask_init(iim_out, jjm_out, nbland_out)
    !
    ! ARGUMENTS
    !
    INTEGER(i_std), INTENT(out) :: iim_out, jjm_out, nbland_out
    !
    ! LOCAL
    !
    CHARACTER(LEN=120)                                  :: filename                !! File name for topography (var: ROSE in "etopo20.nc")
    CHARACTER(LEN=120)                                  :: contfracfile            !! filename for contfrac
    INTEGER(i_std)                                      :: iim_full, jjm_full      !! Lon/Lat dimensions for ROSE in etopo20.nc 
    REAL(r_std), ALLOCATABLE, DIMENSION (:,:)           :: lon_full, lat_full      !! Lon/Lat array for ROSE
    REAL(r_std), ALLOCATABLE, DIMENSION (:,:)           :: etopo, contfrac_full    !! Array for reading ROSE in etopo20.nc
    INTEGER(i_std), DIMENSION(1)                        :: ibeg, iend
    INTEGER(i_std), DIMENSION(1)                        :: jbeg, jend
    INTEGER(i_std)                                      :: fid, i, j
    REAL(r_std)                                         :: dx, dy
    !
    REAL(r_std), DIMENSION(2)                           :: zoom_lon, zoom_lat               !! For crop region when reading
    !
    LOGICAL                                             :: buildcontfrac=.FALSE.
    !
    !  Read the full topography which we will use to get the land-see mask for all resolutions.
    !
    !Config Key   = TOPOFILE
    !Config Desc  = File containing high resolution topographic informations
    !Config If    = 
    !Config Def   = etopo20.nc
    !Config Help  = This needs to be a netCDF file.
    !Config Units = [-] 
    !- 
    filename = 'NONE'
    CALL getin('TOPOFILE', filename)
    WRITE(*,*) "TOPOFILE = ", filename
    !
    !Config Key   = CONTFRACFILE
    !Config Desc  = File containing a contfrac variable as defined in ORCHIDEE.
    !Config If    = 
    !Config Def   = NONE
    !Config Help  = This needs to be a netCDF file and can be any history file ORCHIDEE
    !               has produced ina previous simulation.
    !Config Units = [-] 
    !- 
    contfracfile = 'NONE'
    CALL getin('CONTFRACFILE', contfracfile)
    WRITE(*,*) "CONTFRACFILE = ", contfracfile
    !
    ! Initialize Orchidee parameter and calendar
    !
    IF ( INDEX(filename,"NONE") <= 0 ) THEN
       IF ( is_root_prc) THEN
          CALL topo_getsize(filename, iim_full, jjm_full)
          !
          ALLOCATE(lon_full(iim_full,jjm_full),lat_full(iim_full,jjm_full))
          ALLOCATE(etopo(iim_full,jjm_full))
          !
          CALL topo_getvar(filename, iim_full, jjm_full, lon_full, lat_full, etopo)
          !
          WRITE(*,*) MINVAL(lon_full), ' < LON_FULL < ', MAXVAL(lon_full)
          WRITE(*,*) MINVAL(lat_full), ' < LAT_FULL < ', MAXVAL(lat_full)
          WRITE(*,*) 'ETOPO :', MINVAL(etopo), MAXVAL(etopo)
          !
          buildcontfrac=.TRUE.
          !
       ENDIF
    ELSE IF ( INDEX( contfracfile,"NONE") <= 0 ) THEN
       IF ( is_root_prc) THEN
          !
          ! Read lon, lat and contfrac from a netCDF file generated by ORCHIDEE
          !
          CALL opencontfrac(contfracfile, iim_full, jjm_full, fid)
          !
          ALLOCATE(lon_full(iim_full,jjm_full), lat_full(iim_full,jjm_full), contfrac_full(iim_full,jjm_full))
          CALL readcontfrac(fid, contfracfile, iim_full, jjm_full, lon_full, lat_full, contfrac_full)
          !
          buildcontfrac=.FALSE.
          !
       ENDIF
    ELSE
       WRITE(*,*) "Neither a orography nore a contfrac file was provided"
       WRITE(*,*) "CONTFRACFILE = ", contfracfile
       write(*,*) "TOPOFILE = ", filename
    ENDIF
    !
    !-
    !- Define the zoom
    !-
    zoom_lon=(/-180,180/)
    zoom_lat=(/-90,90/)
    !
    !Config Key   = LIMIT_WEST
    !Config Desc  = Western limit of region
    !Config If    = [-]
    !Config Def   = -180.
    !Config Help  = Western limit of the region we are 
    !Config         interested in. Between -180 and +180 degrees
    !Config         The model will use the smalest regions from
    !Config         region specified here and the one of the forcing file.
    !Config Units = [Degrees] 
    !- 
    CALL getin_p('LIMIT_WEST',zoom_lon(1))
    !-
    !Config Key   = LIMIT_EAST
    !Config Desc  = Eastern limit of region
    !Config If    = [-]
    !Config Def   = 180.
    !Config Help  = Eastern limit of the region we are
    !Config         interested in. Between -180 and +180 degrees
    !Config         The model will use the smalest regions from
    !Config         region specified here and the one of the forcing file.
    !Config Units = [Degrees] 
    !-
    CALL getin_p('LIMIT_EAST',zoom_lon(2))
    !-
    !Config Key   = LIMIT_NORTH
    !Config Desc  = Northern limit of region
    !Config If    = [-]
    !Config Def   = 90.
    !Config Help  = Northern limit of the region we are
    !Config         interested in. Between +90 and -90 degrees
    !Config         The model will use the smalest regions from
    !Config         region specified here and the one of the forcing file.
    !Config Units = [Degrees]
    !-
    CALL getin_p('LIMIT_NORTH',zoom_lat(2))
    !-
    !Config Key   = LIMIT_SOUTH
    !Config Desc  = Southern limit of region
    !Config If    = [-]
    !Config Def   = -90.
    !Config Help  = Southern limit of the region we are
    !Config         interested in. Between 90 and -90 degrees
    !Config         The model will use the smalest regions from
    !Config         region specified here and the one of the forcing file.
    !Config Units = [Degrees]
    !-
    CALL getin_p('LIMIT_SOUTH',zoom_lat(1))
    IF ( (zoom_lon(1)+180 < EPSILON(zoom_lon(1))) .AND. (zoom_lon(2)-180 < EPSILON(zoom_lon(2))) .AND.&
        &(zoom_lat(1)+90 < EPSILON(zoom_lat(1))) .AND. (zoom_lat(2)-90 < EPSILON(zoom_lat(2))) ) THEN
       !
       !Config Key   = WEST_EAST
       !Config Desc  = Longitude interval to use from the forcing data
       !Config If    = [-]
       !Config Def   = -180, 180
       !Config Help  = This function allows to zoom into the forcing data
       !Config Units = [degrees east] 
       !- 
       CALL getin('WEST_EAST', zoom_lon)
       !
       !Config Key   = SOUTH_NORTH
       !Config Desc  = Latitude interval to use from the forcing data
       !Config If    = [-]
       !Config Def   = -90, 90
       !Config Help  = This function allows to zoom into the forcing data
       !Config Units = [degrees north] 
       !- 
       CALL getin('SOUTH_NORTH', zoom_lat)
    ENDIF
    !
    ! Get the finer coarser grid
    !
    !
    ! We see if the user has specified the resolution, in number of
    ! points he wishes to use.
    !
    IF ( buildcontfrac ) THEN
       IF ( is_root_prc) THEN
          iim=-1
          jjm=-1
          !Config Key   = IIM_LON
          !Config Desc  = Number of points in longitude
          !Config If    = [-]
          !Config Def   = -1
          !Config Help  = This allows to select the resolution at which the land/sea mask should be built.
          !Config Units = - 
          CALL getin('IIM_LON', iim)
          !Config Key   = IIM_LAT
          !Config Desc  = Number of points in latitude
          !Config If    = [-]
          !Config Def   = -1
          !Config Help  = This allows to select the resolution at which the land/sea mask should be built.
          !Config Units = - 
          CALL getin('JJM_LAT', jjm)
          !
          IF (iim > 0 .AND. jjm > 0 ) THEN
             !
             ALLOCATE(lon(iim,jjm), lat(iim,jjm), orog(iim,jjm))
             ! Trung: Calculate resolution
             dx = (MAXVAL(zoom_lon)-MINVAL(zoom_lon))/REAL(iim, r_std)
             dy = (MAXVAL(zoom_lat)-MINVAL(zoom_lat))/REAL(jjm, r_std)
             WRITE (*,*) "Resolution in longitude [degrees] dx = ",dx
             WRITE (*,*) "Resolution in laltiude [degrees] dy = ",dy
             !
             ! Generate lon and lat
             !
             DO i=1,iim
                lon(i,:) = MINVAL(zoom_lon) + (i-1)*dx + dx/2.0 
             ENDDO
             DO j=1,jjm
                lat(:,j) = MINVAL(zoom_lat) + (j-1)*dy + dy/2.0
             ENDDO
             !
             ! Interpolate to new orography
             !
             CALL interpol(iim_full, jjm_full, lon_full, lat_full, etopo, iim, jjm, lon, lat, orog)
             !
             ! Write orog to netcdf for checking (deleted)
             !
          ELSE
             ! Just transfer the data over the zoomed area
             ibeg=MINLOC(ABS(lon_full(:,1)-MINVAL(zoom_lon)))
             iend=MINLOC(ABS(lon_full(:,1)-MAXVAL(zoom_lon)))
             jbeg=MINLOC(ABS(lat_full(1,:)-MINVAL(zoom_lat)))
             jend=MINLOC(ABS(lat_full(1,:)-MAXVAL(zoom_lat)))
             iim = (iend(1)-ibeg(1))+1
             jjm = (jend(1)-jbeg(1))+1
             !
             !
             ALLOCATE(lon(iim,jjm), lat(iim,jjm), orog(iim,jjm), contfrac(iim,jjm))
             !
             DO j=1,jjm
                lon(:,j) = lon_full(ibeg(1):iend(1),j)
             ENDDO
             DO i=1,iim
                lat(i,:) = lat_full(i,jbeg(1):jend(1))
             ENDDO
             ! Some treatment for the orogography
             orog(:,:) = 0.0
             DO i=ibeg(1),iend(1)
                DO j=jbeg(1),jend(1)
                   orog(i-ibeg(1)+1,j-jbeg(1)+1) = MAX(etopo(i,j), 0.0)
                ENDDO
             ENDDO
             !
             ! Generate the contfrac field
             !
             contfrac(:,:)= 0.0
             DO i=1,iim
                DO j=1,jjm
                   IF ( orog(i,j) > 0 ) THEN
                      contfrac(i,j) = 1.0
                   ENDIF
                ENDDO
             ENDDO
          ENDIF
          !
       ENDIF
       !
    ELSE
       !
       ! Just zoom into the contfrac if needed and build a dummy orog
       !
       IF ( is_root_prc) THEN
          IF ( MINVAL(zoom_lon) > -180 .OR. MAXVAL(zoom_lon) < 180 .OR.&
               & MINVAL(zoom_lat) > -90 .OR. MAXVAL(zoom_lat) < 90 ) THEN
             ! Just transfer the data over the zoomed area
             ibeg=MINLOC(ABS(lon_full(:,1)-MINVAL(zoom_lon)))
             iend=MINLOC(ABS(lon_full(:,1)-MAXVAL(zoom_lon)))
             jbeg=MINLOC(ABS(lat_full(1,:)-MINVAL(zoom_lat)))
             jend=MINLOC(ABS(lat_full(1,:)-MAXVAL(zoom_lat)))
             iim = (iend(1)-ibeg(1))+1
             jjm = (jend(1)-jbeg(1))+1
             !
             !
             ALLOCATE(lon(iim,jjm), lat(iim,jjm), orog(iim,jjm), contfrac(iim,jjm))
             !
             DO j=1,jjm
                lon(:,j) = lon_full(ibeg(1):iend(1),j)
             ENDDO
             DO i=1,iim
                lat(i,:) = lat_full(i,jbeg(1):jend(1))
             ENDDO
             ! Some treatment for the orogography
             orog(:,:) = 0.0
             contfrac(:,:) = 0.0
             DO i=ibeg(1),iend(1)
                DO j=jbeg(1),jend(1)
                   orog(i-ibeg(1)+1,j-jbeg(1)+1) = 1.0
                   contfrac(i-ibeg(1)+1,j-jbeg(1)+1) = contfrac_full(i,j)
                ENDDO
             ENDDO
             !
          ELSE
             iim = iim_full
             jjm = jjm_full
             ALLOCATE(lon(iim,jjm), lat(iim,jjm), orog(iim,jjm), contfrac(iim,jjm))
             lon(:,:) = lon_full(:,:)
             lat(:,:) = lat_full(:,:)
             contfrac(:,:) = contfrac_full(:,:)
             orog(:,:) = contfrac_full(:,:)+100
          ENDIF
       ENDIF
       !
    ENDIF
    !
    IF ( is_root_prc) THEN
       WRITE(*,*) "Dimensions : ", iim, jjm
       WRITE(*,*) MINVAL(lon), ' < LON < ', MAXVAL(lon)
       WRITE(*,*) MINVAL(lat), ' < LAT < ', MAXVAL(lat)  
       WRITE(*,*) MINVAL(orog), " < OROG  < ", MAXVAL(orog)
       WRITE(*,*) MINVAL(contfrac), " < contfrac  < ", MAXVAL(contfrac)
    ENDIF
    !
    ! Free some memory
    !
    IF (ALLOCATED(lon_full)) DEALLOCATE(lon_full)
    IF (ALLOCATED(lat_full)) DEALLOCATE(lat_full)
    IF (ALLOCATED(etopo)) DEALLOCATE(etopo)
    IF (ALLOCATED(contfrac_full)) DEALLOCATE(contfrac_full)
    !
    ! Distribute the information to all processors
    !
    CALL bcast(iim)
    CALL bcast(jjm)
    IF ( .NOT. ALLOCATED(lon)) ALLOCATE(lon(iim,jjm))
    IF ( .NOT. ALLOCATED(lat)) ALLOCATE(lat(iim,jjm))
    IF ( .NOT. ALLOCATED(orog)) ALLOCATE(orog(iim,jjm))
    IF ( .NOT. ALLOCATED(contfrac)) ALLOCATE(contfrac(iim,jjm))
    CALL bcast(lon)
    CALL bcast(lat)
    CALL bcast(orog)
    CALL bcast(contfrac)
    !
    ! Count the number of land points
    !
    nbland = COUNT(contfrac > 0.0)
    !
    iim_out = iim
    jjm_out = jjm
    nbland_out = nbland
    !
  END SUBROUTINE getlandseamask_init
  !!
  !================================================================================================================================
  !! SUBROUTINE   : getlandseamask_read
  !!
  !> BRIEF         Returns the land/sea mask
  !!
  !! DESCRIPTION  : None
  !!
  !! RECENT CHANGE(S): None
  !!
  !! MAIN OUTPUT VARIABLE(S):
  !!
  !! REFERENCES   : None
  !!
  !! FLOWCHART    : None
  !! \n
  !_================================================================================================================================
  !!
  SUBROUTINE getlandseamask_read(lon_out, lat_out, mask_out, orog_out)
    !
    ! ARGUMENTS
    !
    REAL(r_std), DIMENSION(:,:), INTENT(out) :: lon_out, lat_out
    REAL(r_std), DIMENSION(:,:), INTENT(out) :: mask_out, orog_out
    !
    ! LOCAL
    !
    lon_out(:,:) = lon(:,:)
    lat_out(:,:) = lat(:,:)
    mask_out(:,:) = contfrac(:,:)
    orog_out(:,:) = orog(:,:)
    !
  END SUBROUTINE getlandseamask_read
  !!
  !!
  !
  !================================================================================================================================
  !! SUBROUTINE   : opencontfrac
  !!
  !> BRIEF         This subroutine gets the dimensions of the fields in the filename. Lat and lon in particular.
  !!
  !! DESCRIPTION  : This routines caters for ORCHIDEE as well as the WRF geo files.
  !!
  !! RECENT CHANGE(S): None
  !!
  !! MAIN OUTPUT VARIABLE(S):
  !!
  !! REFERENCES   : None
  !!
  !! FLOWCHART    : None
  !! \n
  !_================================================================================================================================
  SUBROUTINE opencontfrac(filename, iim_full, jjm_full, fid)
    !
    ! ARGUMENTS
    !
    CHARACTER(LEN=*), INTENT(in) :: filename
    INTEGER(i_std), INTENT(out) :: iim_full, jjm_full, fid
    !
    INTEGER(i_std) :: iret, i, lll, ndims, nvars
    CHARACTER(LEN=20) :: dimname
    !
    iret = NF90_OPEN(filename, NF90_NOWRITE, fid)
    IF (iret /= NF90_NOERR) THEN
       WRITE(*,*) "==>",trim(nf90_strerror(iret))
       WRITE(*,*) "Error opening ", filename
       STOP
    ENDIF
    iret = NF90_INQUIRE (fid, nDimensions=ndims, nVariables=nvars)
    DO i=1,ndims
       iret = NF90_INQUIRE_DIMENSION (fid, i, name=dimname, len=lll)
       CALL change_to_lowercase(dimname)
       IF ( INDEX(dimname,"lon") > 0 ) THEN
          iim_full = lll
       ELSE IF (INDEX(dimname,"west_east") > 0 .AND. LEN_TRIM(dimname) == LEN_TRIM("west_east")) THEN
          iim_full = lll
       ELSE IF ( INDEX(dimname,"lat") > 0 ) THEN
          jjm_full = lll
       ELSE IF (INDEX(dimname,"south_north") > 0 .AND. LEN_TRIM(dimname) == LEN_TRIM("south_north")) THEN
          jjm_full = lll
       ENDIF
    ENDDO
    !
    WRITE(*,*) "XXX iimf, jjmf = ", iim_full, jjm_full
    !
  END SUBROUTINE opencontfrac

  !================================================================================================================================
  !! SUBROUTINE   : readcontfrac
  !!
  !> BRIEF         This subroutine gets the topography out of the file
  !!
  !! DESCRIPTION  : This routines caters for ORCHIDEE as well as the WRF geo files.
  !!
  !! RECENT CHANGE(S): None
  !!
  !! MAIN OUTPUT VARIABLE(S):
  !!
  !! REFERENCES   : None
  !!
  !! FLOWCHART    : None
  !! \n
  !_================================================================================================================================
  SUBROUTINE readcontfrac(fid, filename, iimf, jjmf, lon, lat, contf)
    !
    ! ARGUMENTS
    !
    INTEGER(i_std), INTENT(in)               :: fid
    CHARACTER(LEN=*), INTENT(in)             :: filename
    INTEGER(i_std), INTENT(in)               :: iimf, jjmf
    REAL(r_std), DIMENSION(:,:), INTENT(out) :: lon, lat
    REAL(r_std), DIMENSION(:,:), INTENT(out) :: contf
    !
    ! LOCAL
    !
    INTEGER(i_std) :: iret, ndims, nvars, xid, yid, vid, wid
    INTEGER(i_std) :: xndims, yndims, vndims, wndims
    INTEGER(i_std) :: iv, ndimsvar, i, j
    INTEGER(i_std), DIMENSION(4) :: dimids
    CHARACTER(LEN=60) :: varname, units
    CHARACTER(LEN=3)  :: model='ORC'
    INTEGER(i_std), DIMENSION(3) :: start, count
    REAL(r_std), ALLOCATABLE, DIMENSION(:) :: vtmp
    !
    iret = NF90_INQUIRE (fid, nDimensions=ndims, nVariables=nvars)
    !
    xid = -1
    yid = -1
    vid = -1
    wid = -1
    !
    DO iv=1,nvars
       !
       iret = NF90_INQUIRE_VARIABLE(fid, iv, name=varname, ndims=ndimsvar, dimids=dimids)
       iret = NF90_GET_ATT(fid, iv, 'units', units)
       CALL change_to_lowercase(units)
       CALL change_to_lowercase(varname)
       !
       IF ( INDEX(units,"east") > 0 .OR. INDEX(varname, "xlong_m") > 0 ) THEN
          xid = iv
          xndims = ndimsvar
          IF ( INDEX(varname, "xlong_m") > 0 ) THEN
             model='WRF'
          ENDIF
       ELSE IF (INDEX(units,"north") > 0 .OR. INDEX(varname, "xlat_m") > 0 ) THEN
          yid = iv
          yndims = ndimsvar
          IF ( INDEX(varname, "xlat_m") > 0 ) THEN
             model='WRF'
          ENDIF
       ENDIF
       !
       IF ( INDEX(varname,"contfrac") > 0 ) THEN
          vid = iv
          vndims = ndimsvar
       ENDIF
       IF ( INDEX(varname,"landmask") > 0 ) THEN
          wid = iv
          wndims = ndimsvar
       ENDIF
       !
    ENDDO
    !
    ! Get longitude
    !
    IF ( xid > 0 ) THEN
       IF ( xndims == 1 ) THEN
          ALLOCATE(vtmp(iimf))
          iret = NF90_GET_VAR(fid, xid, vtmp)
          DO j=1,jjmf
             lon(:,j) = vtmp(:)
          ENDDO
          DEALLOCATE(vtmp)
       ELSE IF ( xndims == 2 ) THEN
          iret = NF90_GET_VAR(fid, xid, lon)
       ELSE IF ( xndims == 3 ) THEN
          start = (/1,1,1/)
          count = (/iimf,jjmf,1/)
          iret = NF90_GET_VAR(fid, xid, lon, start, count)
       ELSE
          WRITE(*,*) "Unforeseen number of dimensions for longitude ", xndims
          STOP
       ENDIF
    ELSE
       WRITE(*,*) "could not find the longitude in file ", filename
       STOP
    ENDIF
    !
    ! Get latitude
    !
    IF ( yid > 0 ) THEN
       IF ( yndims == 1 ) THEN
          ALLOCATE(vtmp(jjmf))
          iret = NF90_GET_VAR(fid, xid, vtmp)
          DO i=1,iimf
             lat(i,:) = vtmp(:)
          ENDDO
          DEALLOCATE(vtmp)
       ELSE IF ( yndims == 2 ) THEN
          iret = NF90_GET_VAR(fid, yid, lat)
       ELSE IF ( yndims == 3 ) THEN
          start = (/1,1,1/)
          count = (/iimf,jjmf,1/)
          iret = NF90_GET_VAR(fid, yid, lat, start, count)
       ELSE
          WRITE(*,*) "Unforeseen number of dimensions for latitude ", yndims
          STOP
       ENDIF
    ELSE
       WRITE(*,*) "could not find the latitude in file ", filename
       STOP
    ENDIF
    !
    ! Get the variable
    !
    contf(:,:) = 0.0
    IF ( model == "ORC" ) THEN
       IF ( vid > 0 ) THEN
          iret = NF90_GET_VAR(fid, vid, contf)
          IF (iret /= NF90_NOERR) THEN
             WRITE(*,*) "==>",trim(nf90_strerror(iret))
             WRITE(*,*) "Cannot read variable contfrac"
             STOP
          ENDIF
       ELSE
          WRITE(*,*) "could not find the variable contfrac in file ", filename
          STOP
       ENDIF
    ELSE IF ( model == "WRF" ) THEN
       !
       ! This is a time varying variable so start and count needs to be specified
       !
       start = (/1,1,1/)
       count = (/iimf,jjmf,1/)
       IF ( wid > 0 ) THEN
          iret = NF90_GET_VAR(fid, wid, contf, start, count)
          IF (iret /= NF90_NOERR) THEN
             WRITE(*,*) "==>",trim(nf90_strerror(iret))
             WRITE(*,*) "Cannot read variable landmask"
             STOP
          ENDIF
       ELSE
          WRITE(*,*) "could not find the variable contfrac in file ", filename
          STOP
       ENDIF
       WRITE(*,*) "XXX", MINVAL(contf), " < contf < ", MAXVAL(contf)
    ELSE
       WRITE(*,*) "Unknown model = ", model
       STOP
    ENDIF
    !
    ! Delete undef values
    !
    DO i=1,iimf
       DO j=1,jjmf
          IF ( contf(i,j) > 1.0 .OR. contf(i,j) < 0.0 ) THEN
             contf(i,j) = 0.0
          ENDIF
       ENDDO
    ENDDO
    !
    iret = NF90_CLOSE(fid)
    !
  END SUBROUTINE readcontfrac
  !!
  !
  !================================================================================================================================
  !! SUBROUTINE   : topo_getsize
  !!
  !> BRIEF         This subroutine gets the size of the topography field
  !!
  !! DESCRIPTION  : None
  !!
  !! RECENT CHANGE(S): None
  !!
  !! MAIN OUTPUT VARIABLE(S):
  !!
  !! REFERENCES   : None
  !!
  !! FLOWCHART    : None
  !! \n
  !_================================================================================================================================
  SUBROUTINE topo_getsize(filename, iim_full, jjm_full)
    !
    USE defprec
    USE netcdf
    !
    ! ARGUMENTS
    !
    CHARACTER(LEN=*), INTENT(in) :: filename
    INTEGER(i_std), INTENT(out)  :: iim_full, jjm_full
    !
    ! LOCAL
    !
    INTEGER(i_std) :: fid, iret, xid, yid, ndims, nvars, i, lll
    CHARACTER(LEN=20) :: dimname
    !
    iret = NF90_OPEN(filename, NF90_NOWRITE, fid)
    IF (iret /= NF90_NOERR) THEN
       WRITE(*,*) "Error opening ", filename
       STOP
    ENDIF
    iret = NF90_INQUIRE (fid, nDimensions=ndims, nVariables=nvars)
    DO i=1,ndims
       iret = NF90_INQUIRE_DIMENSION (fid, i, name=dimname, len=lll)
       CALL change_to_lowercase(dimname)
       IF ( INDEX(dimname,"x") > 0 ) THEN
          iim_full = lll
       ELSE IF ( INDEX(dimname,"y") > 0 ) THEN
          jjm_full = lll
       ENDIF
    ENDDO
    !
    iret = NF90_CLOSE(fid)
    !
  END SUBROUTINE topo_getsize
  !================================================================================================================================
  !! SUBROUTINE   : topo_getvar
  !!
  !> BRIEF         This subroutine gets the topography out of the file
  !!
  !! DESCRIPTION  : None
  !!
  !! RECENT CHANGE(S): None
  !!
  !! MAIN OUTPUT VARIABLE(S):
  !!
  !! REFERENCES   : None
  !!
  !! FLOWCHART    : None
  !! \n
  !_================================================================================================================================
  SUBROUTINE topo_getvar(filename, iim, jjm, lon, lat, etopo)
    !
    USE defprec
    USE netcdf
    !
    ! ARGUMENTS
    !
    CHARACTER(LEN=*), INTENT(in)                 :: filename
    INTEGER(i_std), INTENT(in)                   :: iim, jjm
    REAL(r_std), DIMENSION(iim,jjm), INTENT(out) :: lon 
    REAL(r_std), DIMENSION(iim,jjm), INTENT(out) :: lat
    REAL(r_std), DIMENSION(iim,jjm), INTENT(out) :: etopo
    !
    ! LOCAL
    !
    INTEGER(i_std) :: fid, iret, xid, yid, vid
    INTEGER(i_std) :: i, j, iv, ndimsvar, ndims, nvars
    CHARACTER(LEN=20) :: varname, units
    INTEGER(i_std), DIMENSION(1) :: il
    REAL(r_std), ALLOCATABLE, DIMENSION(:)   :: lon_loc, lon_read, lat_read
    REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: etopo_loc
    INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: dimids, londim_ids, latdim_ids
    !
    ALLOCATE(lon_loc(iim))
    ALLOCATE(lon_read(iim))
    ALLOCATE(lat_read(jjm))
    ALLOCATE(etopo_loc(iim,jjm))
    ALLOCATE(dimids(NF90_MAX_VAR_DIMS), londim_ids(NF90_MAX_VAR_DIMS), latdim_ids(NF90_MAX_VAR_DIMS))
    !
    iret = NF90_OPEN(filename, NF90_NOWRITE, fid)
    IF (iret /= NF90_NOERR) THEN
       WRITE(*,*) "Error opening ", filename
       STOP
    ENDIF
    !
    iret = NF90_INQUIRE (fid, nDimensions=ndims, nVariables=nvars)
    !
    ! Identify the variables
    !
    xid = -1
    yid = -1
    vid = -1
    !
    DO iv=1,nvars
       iret = NF90_INQUIRE_VARIABLE(fid, iv, name=varname, ndims=ndimsvar, dimids=dimids)
       iret = NF90_GET_ATT(fid, iv, 'units', units)
       CALL change_to_lowercase(units)
       !
       IF ( INDEX(units,"east") > 0 ) THEN
          xid = iv
       ELSE IF (INDEX(units,"north") > 0 ) THEN
          yid = iv
       ELSE IF ( INDEX(units,"meters") > 0 ) THEN
          vid = iv
       ENDIF
    ENDDO
    !
    ! Get variables now that they have been identified
    !
    IF ( xid > 0 ) THEN
       iret = NF90_GET_VAR(fid, xid, lon_read)
    ELSE
       WRITE(*,*) "cound not find the longitude in file ", filename
       STOP
    ENDIF
    IF ( yid > 0 ) THEN
       iret = NF90_GET_VAR(fid, yid, lat_read)
    ELSE
       WRITE(*,*) "cound not find the latitude in file ", filename
       STOP
    ENDIF
    IF ( vid > 0 ) THEN
       iret = NF90_GET_VAR(fid, vid, etopo_loc)
    ELSE
       WRITE(*,*) "cound not find the topography in file ", filename
       STOP
    ENDIF
    iret = NF90_CLOSE(fid)
    !
    ! Shift the topography around the longitude to get a -180:180 range 
    ! Only if needed !
    !
    IF ( MAXVAL(lon_read) > 180. ) THEN
       ! 
       lon_loc(:) = lon_read(:)
       !
       DO i=1,iim
          IF (lon_loc(i) > 180. ) THEN
             lon_loc(i) = lon_loc(i)-360.
          ENDIF
       ENDDO
       !
       DO i=1,iim
          il = MINLOC(lon_loc)
          lon_read(i) = lon_loc(il(1))
          etopo(i,:) = etopo_loc(il(1),:)
          lon_loc(il(1)) = 9999999.999999
       ENDDO
       !
       DO j=1,jjm
          lon(:,j) = lon_read(:)
       ENDDO
       DO i=1,iim
          lat(i,:) = lat_read(:)
       ENDDO
    ELSE
       DO j=1,jjm
          lon(:,j) = lon_read(:)
       ENDDO
       DO i=1,iim
          lat(i,:) = lat_read(:)
       ENDDO
       etopo(:,:) = etopo_loc(:,:)
    ENDIF
    !
    !
    DEALLOCATE(lon_read)
    DEALLOCATE(lon_loc)
    DEALLOCATE(lat_read)
    DEALLOCATE(etopo_loc)
    !
  END SUBROUTINE topo_getvar
  !
  !================================================================================================================================
  !! SUBROUTINE   : interpol
  !!
  !> BRIEF         This subroutine interpolates topography to new resolution.
  !!
  !! DESCRIPTION  : None
  !!
  !! RECENT CHANGE(S): None
  !!
  !! MAIN OUTPUT VARIABLE(S):
  !!
  !! REFERENCES   : None
  !!
  !! FLOWCHART    : None
  !! \n
  !_================================================================================================================================

  SUBROUTINE interpol(iim_fine, jjm_fine, lon_fine, lat_fine, orog, iim, jjm, lon, lat, neworog)
    !
    USE constantes_var
    USE ioipsl_para
    !
    IMPLICIT NONE
    ! Parameter
    !INTEGER(i_std), PARAMETER                                  :: nbvmax = 800
    INTEGER(i_std), PARAMETER                                  :: nbvmax = 30000
    !
    INTEGER(i_std)                                             :: iim_fine, jjm_fine, iim, jjm
    REAL(r_std), DIMENSION (iim_fine,jjm_fine)                 :: lon_fine
    REAL(r_std), DIMENSION (iim_fine,jjm_fine)                 :: lat_fine
    REAL(r_std), DIMENSION (iim_fine,jjm_fine)                 :: orog
    REAL(r_std), DIMENSION (iim,jjm)                           :: lon
    REAL(r_std), DIMENSION (jjm,iim)                           :: lat
    REAL(r_std), DIMENSION (iim,jjm)                           :: neworog
    !
    REAL(r_std), DIMENSION (nbvmax)                            :: area, topo
    !
    INTEGER(i_std)                                             :: ip, jp, ig, jg, i, fopt, lastjp
    REAL(r_std)                                                :: dx, dy, ax, ay, ox, lon_up, lon_low, lat_up, lat_low, sgn
    REAL(r_std)                                                :: totarea, landarea, height
    REAL(r_std), ALLOCATABLE, DIMENSION (:,:)                  :: lon_ful, lat_ful, laup_rel, loup_rel, lalow_rel, lolow_rel
    REAL(r_std)                                                :: lonrel, lolowrel, louprel, coslat
    !
    ! Allocate arrays
    !
    ALLOCATE (laup_rel(iim_fine,jjm_fine))
    ALLOCATE (loup_rel(iim_fine,jjm_fine))
    ALLOCATE (lalow_rel(iim_fine,jjm_fine))
    ALLOCATE (lolow_rel(iim_fine,jjm_fine))
    ALLOCATE (lat_ful(iim_fine+2,jjm_fine+2))
    ALLOCATE (lon_ful(iim_fine+2,jjm_fine+2))
    !
    ! Duplicate the border assuming we have a global grid going from west to east
    !
    DO jp=2,jjm_fine+1
       lon_ful(2:iim_fine+1,jp) = lon_fine(1:iim_fine,jp)
    ENDDO
    DO ip=2,iim_fine+1
       lat_ful(ip,2:jjm_fine+1) = lat_fine(ip,1:jjm_fine)
    ENDDO
    !
    IF ( lon_fine(iim_fine,1) .LT. lon_ful(2,2)) THEN
       lon_ful(1,2:jjm_fine+1) = lon_fine(iim_fine,:)
       lat_ful(1,2:jjm_fine+1) = lat_fine(1,1:jjm_fine)
    ELSE
       lon_ful(1,2:jjm_fine+1) = lon_fine(iim_fine,:)-360.
       lat_ful(1,2:jjm_fine+1) = lat_fine(1,1:jjm_fine)
    ENDIF

    IF ( lon_fine(1,1) .GT. lon_ful(iim_fine+1,2)) THEN
       lon_ful(iim_fine+2,2:jjm_fine+1) = lon_fine(1,1:jjm_fine)
       lat_ful(iim_fine+2,2:jjm_fine+1) = lat_fine(iim_fine,1:jjm_fine)
    ELSE
       lon_ful(iim_fine+2,2:jjm_fine+1) = lon_fine(1,1:jjm_fine)+360
       lat_ful(iim_fine+2,2:jjm_fine+1) = lat_fine(iim_fine,1:jjm_fine)
    ENDIF
    !
    sgn = lat_fine(1,1)/ABS(lat_fine(1,1))
    lat_ful(2:iim_fine+1,1) = sgn*180 - lat_fine(1,1)
    sgn = lat_fine(iim_fine,jjm_fine)/ABS(lat_fine(iim_fine,jjm_fine))
    lat_ful(2:iim_fine+1,jjm_fine+2) = sgn*180 - lat_fine(iim_fine,jjm_fine)
    lat_ful(1,1) = lat_ful(iim_fine+1,1)
    lat_ful(iim_fine+2,1) = lat_ful(2,1)
    lat_ful(1,jjm_fine+2) = lat_ful(iim_fine+1,jjm_fine+2)
    lat_ful(iim_fine+2,jjm_fine+2) = lat_ful(2,jjm_fine+2)
    !
    ! Add the longitude lines to the top and bottom
    !
    lon_ful(:,1) = lon_ful(:,2) 
    lon_ful(:,jjm_fine+2) = lon_ful(:,jjm_fine+1) 
    !
    ! Get the upper and lower limits of each grid box
    !
    DO ip=1,iim_fine
       DO jp=1,jjm_fine
          loup_rel(ip,jp) =MAX(0.5*(lon_ful(ip,jp+1)+lon_ful(ip+1,jp+1)), 0.5*(lon_ful(ip+1,jp+1)+lon_ful(ip+2,jp+1)))
          lolow_rel(ip,jp) =MIN(0.5*(lon_ful(ip,jp+1)+lon_ful(ip+1,jp+1)), 0.5*(lon_ful(ip+1,jp+1)+lon_ful(ip+2,jp+1)))
          laup_rel(ip,jp) =MAX(0.5*(lat_ful(ip+1,jp)+lat_ful(ip+1,jp+1)), 0.5*(lat_ful(ip+1,jp+1)+lat_ful(ip+1,jp+2)))
          lalow_rel(ip,jp) =MIN(0.5*(lat_ful(ip+1,jp)+lat_ful(ip+1,jp+1)), 0.5*(lat_ful(ip+1,jp+1)+lat_ful(ip+1,jp+2)))
       ENDDO
    ENDDO
    !
    ! Now we take each grid point and find out which values from the forcing we need to average
    !
    dx = ABS(lon(2,1) - lon(1,1))
    dy = ABS(lat(1,2) - lat(1,1))
    !
    DO ig = 1,iim
       DO jg = 1,jjm
          !
          !
          lon_up = lon(ig,jg) + dx/2.0
          lon_low = lon(ig,jg) - dx/2.0
          !
          !
          lat_up = lat(ig,jg) + dy/2.0
          lat_low = lat(ig,jg) - dy/2.0
          !
          !  Find the grid boxes from the data that go into the model's boxes
          !  We still work as if we had a regular grid ! Well it needs to be localy regular so
          !  so that the longitude at the latitude of the last found point is close to the one of the next point.
          !
          totarea = 0.0
          landarea = 0.0
          height = 0.0
          !
          fopt = 0
          lastjp = 1
          DO ip=1,iim_fine
             !
             DO jp = 1, jjm_fine
                !  Either the center of the data grid point is in the interval of the model grid or
                !  the East and West limits of the data grid point are on either sides of the border of
                !  the data grid.
                !
                !  To do that correctly we have to check if the grid box sits on the date-line.
                !
                IF ( lon_low < -180.0 ) THEN
                   lonrel = MOD( lon_fine(ip,jp) - 360.0, 360.0)
                   lolowrel = MOD( lolow_rel(ip,jp) - 360.0, 360.0)
                   louprel = MOD( loup_rel(ip,jp) - 360.0, 360.0)
                   !
                ELSE IF ( lon_up > 180.0 ) THEN
                   lonrel = MOD( 360. - lon_fine(ip,jp), 360.0)
                   lolowrel = MOD( 360. - lolow_rel(ip,jp), 360.0)
                   louprel = MOD( 360. - loup_rel(ip,jp), 360.0)
                ELSE
                   lonrel = lon_fine(ip,jp)
                   lolowrel = lolow_rel(ip,jp)
                   louprel = loup_rel(ip,jp)
                ENDIF
                !
                !
                IF ( lonrel > lon_low .AND. lonrel < lon_up .OR. &
                     & lolowrel < lon_low .AND.  louprel > lon_low .OR. &
                     & lolowrel < lon_up  .AND.  louprel > lon_up ) THEN
                   !  
                   !
                   ! Now that we have the longitude let us find the latitude
                   !
                   IF ( lat_fine(ip,jp) > lat_low .AND. lat_fine(ip,jp) < lat_up .OR. &
                        & lalow_rel(ip,jp) < lat_low .AND. laup_rel(ip,jp) > lat_low .OR.&
                        & lalow_rel(ip,jp) < lat_up .AND. laup_rel(ip,jp) > lat_up) THEN
                      !
                      lastjp = jp
                      !
                      fopt = fopt + 1
                      IF ( fopt .GT. nbvmax) THEN
                         WRITE(*,*) 'Please increase nbvmax in subroutine interpol', ig, jg
                         WRITE(*,*) "nbvmax = ", nbvmax
                         STOP
                      ELSE
                         !
                         ! If we sit on the date line we need to do the same transformations as above.
                         ! 
                         IF ( lon_low < -180.0 ) THEN
                            lolowrel = MOD( lolow_rel(ip,jp) - 360.0, 360.0)
                            louprel = MOD( loup_rel(ip,jp) - 360.0, 360.0)
                            !
                         ELSE IF ( lon_up > 180.0 ) THEN
                            lolowrel = MOD( 360. - lolow_rel(ip,jp), 360.0)
                            louprel = MOD( 360. - loup_rel(ip,jp), 360.0)
                         ELSE
                            lolowrel = lolow_rel(ip,jp)
                            louprel = loup_rel(ip,jp)
                         ENDIF
                         !
                         ! Get the area of the fine grid in the model grid
                         !
                         coslat = MAX( COS( lat_fine(ip,jp) * pi/180. ), 0.001 )
                         ax = (MIN(lon_up,louprel)-MAX(lon_low, lolowrel))*pi/180. * R_Earth * coslat
                         ay = (MIN(lat_up, laup_rel(ip,jp))-MAX(lat_low,lalow_rel(ip,jp)))*pi/180. * R_Earth
                         !
                         totarea = totarea + ax*ay
                         IF ( orog(ip, jp) > 0.0 ) THEN
                            landarea = landarea + ax*ay
                            height = height + ax*ay*orog(ip,jp)
                         ENDIF
                         !
                      ENDIF
                   ENDIF
                   !
                ENDIF
                !
             ENDDO
             !
          ENDDO
          !
          IF ( landarea > 0 .AND. landarea > 0.1*totarea ) THEN
             neworog(ig,jg) = height/landarea
          ELSE
             neworog(ig,jg) = 0.0
          ENDIF
          !
       ENDDO
       !
    ENDDO
  END SUBROUTINE interpol

  SUBROUTINE change_to_lowercase (str)
    !---------------------------------------------------------------------
    !- Converts a string into lower case letters
    !---------------------------------------------------------------------
    IMPLICIT NONE
    !-
    CHARACTER(LEN=*) :: str
    !-
    INTEGER :: i,ic
    !---------------------------------------------------------------------
    DO i=1,LEN_TRIM(str)
       ic = IACHAR(str(i:i))
       IF ( (ic >= 65).AND.(ic <= 90) )  str(i:i) = ACHAR(ic+32)
    ENDDO
    !--------------------------
  END SUBROUTINE change_to_lowercase

END MODULE getlandseamask
