grid_atob_m.f90

Go to the documentation of this file.

!*******************************************************************************
!
MODULE grid_atob_m
!
!*******************************************************************************

  USE assert_eq_m, ONLY: assert_eq

  PRIVATE
  PUBLIC :: grille_m, rugosite, sea_ice, rugsoro

CONTAINS

!-------------------------------------------------------------------------------
!
SUBROUTINE fine2coarse(x_i, y_i, x_o, y_o, d_o1, d_i, msk, d_o2)
!
!-------------------------------------------------------------------------------
  IMPLICIT NONE
!-------------------------------------------------------------------------------
! Arguments:
  REAL,              INTENT(IN)  :: x_i(:), y_i(:) !-- INPUT  X&Y COOR. (mi)(ni)
  REAL,              INTENT(IN)  :: x_o(:), y_o(:) !-- OUTPUT X&Y COOR. (mo)(no)
  REAL,              INTENT(OUT) :: d_o1(:,:)      !-- OUTPUT FIELD     (mo,no)
  REAL,    OPTIONAL, INTENT(IN)  :: d_i (:,:)      !-- INPUT FIELD      (mi,ni)
  LOGICAL, OPTIONAL, INTENT(IN)  :: msk (:,:)      !-- MASK             (mo,no)
  REAL,    OPTIONAL, INTENT(OUT) :: d_o2(:,:)      !-- OUTPUT FOR d_i^2 (mo,no)
!-------------------------------------------------------------------------------
! Local variables:
  CHARACTER(LEN=256) :: modname="fine2coarse"
  INTEGER :: mi, ni, ii, ji, mo, no, io, jo, nr(2), m1, n1, m2, n2, mx, my, nn
  LOGICAL :: li, lo
  REAL    :: inc
  INTEGER, DIMENSION(SIZE(x_o),SIZE(y_o)) :: num_tot
  LOGICAL, DIMENSION(SIZE(x_o),SIZE(y_o)) :: found, mask
  REAL,    DIMENSION(SIZE(x_o),SIZE(y_o)) :: dist
  REAL,    DIMENSION(SIZE(x_o))           :: a, b
  REAL,    DIMENSION(SIZE(y_o))           :: c, d
  REAL,    PARAMETER :: thresh=1.e-5
!-------------------------------------------------------------------------------
  mi=SIZE(x_i); ni=SIZE(y_i); mo=SIZE(x_o); no=SIZE(y_o)
  m1=m1; m2=mo; mx=mo; IF(PRESENT(msk)) mx=SIZE(msk,1)
  n1=ni; n2=no; my=no; IF(PRESENT(msk)) my=SIZE(msk,2)
  li=PRESENT(d_i ); IF(li) then; m1=SIZE(d_i ,1); n1=SIZE(d_i ,2); END IF
  lo=PRESENT(d_o2); IF(lo) then; m2=SIZE(d_o2,1); n2=SIZE(d_o2,2); END IF
  mi=assert_eq(mi,m1,trim(modname)//" mi")
  ni=assert_eq(ni,n1,trim(modname)//" ni")
  mo=assert_eq(mo,m2,mx,SIZE(d_o1,1),trim(modname)//" mo")
  no=assert_eq(no,n2,my,SIZE(d_o1,2),trim(modname)//" no")
  mask(:,:)=.true.; IF(PRESENT(msk)) mask(:,:)=msk(:,:)

!--- COMPUTE CELLS INTERFACES COORDINATES OF OUTPUT GRID
  b(mo)=x_o(mo)+(x_o(mo)-x_o(mo-1))/2.0; b(1:mo-1)=(x_o(1:mo-1)+x_o(2:mo))/2.0
  d(no)=y_o(no)+(y_o(no)-y_o(no-1))/2.0; d(1:no-1)=(y_o(1:no-1)+y_o(2:no))/2.0
  a(1 )=x_o(1 )-(x_o(2 )-x_o(1   ))/2.0; a(2:mo  )=   b(1:mo-1)
  c(1 )=y_o(1 )-(y_o(2 )-y_o(1   ))/2.0; c(2:no  )=   d(1:no-1)

!--- ACCUMULATE INPUT POINTS ON OUTPUT GRID
  d_o1(:,:)=0.; num_tot(:,:)=0; inc=1.0
  IF(lo) d_o2(:,:)=0.
  DO ji = 1, ni
    DO ii = 1, mi
      IF(li) inc=d_i(ii,ji)
      DO jo = 1, no
        IF((y_i(ji)-c(jo)<thresh.OR.y_i(ji)-d(jo)>thresh).AND.   &
           (y_i(ji)-c(jo)>thresh.OR.y_i(ji)-d(jo)<thresh)) cycle
        DO io = 1, mo
          IF((x_i(ii)-a(io)<thresh.OR.x_i(ii)-b(io)>thresh).AND. &
             (x_i(ii)-a(io)>thresh.OR.x_i(ii)-b(io)<thresh)) cycle
          num_tot(io,jo)=num_tot(io,jo)+1
          IF(mask(io,jo)) d_o1(io,jo)=d_o1(io,jo)+inc
          IF(.NOT.lo) cycle
          IF(mask(io,jo)) d_o2(io,jo)=d_o2(io,jo)+inc*inc
        END DO
      END DO
    END DO
  END DO

!--- CHECK INPUT POINTS HAVE BEEN FOUND IN EACH OUTPUT CELL
  found(:,:)=num_tot(:,:)/=0
  WHERE(found.AND.mask) d_o1(:,:)=d_o1(:,:)/REAL(num_tot(:,:))
  IF(PRESENT(d_o2)) THEN
    WHERE(found.AND.mask) d_o2(:,:)=d_o2(:,:)/REAL(num_tot(:,:))
    RETURN
  END IF
  nn=count(found); IF(nn==0) RETURN

!--- MISSING POINTS ; USE DISTANCE ON THE SPHERE TO FIND NEAREST POINT nr(2)
  DO io = 1, mo
    DO jo = 1, no
      IF(found(io,jo)) cycle
!      IF(prt_level>=1) PRINT*, "Problem: point out of domain (i,j)=", io,jo
      CALL dist_sphe(x_o(io),y_o(jo),x_i,y_i,dist(:,:))
      nr=minloc(dist(:,:))!; IF(prt_level>=1) PRINT*, "Solution: ", nr
      inc=1.0; IF(li) inc=d_i(nr(1),nr(2))
      IF(mask(io,jo)) d_o1(io,jo)=inc
    END DO
  END DO

END SUBROUTINE fine2coarse
!
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!
SUBROUTINE grille_m(xdata, ydata, entree, x, y, sortie)
!
!-------------------------------------------------------------------------------
! Author:  Z.X. Li (april 1st 1994) (see also A. Harzallah and L. Fairhead)
!-------------------------------------------------------------------------------
! Purpose: Naive method to regrid on a coarser grid.
!   Value at a new point is an average of the old points lcoated in a zone
!   surrounding the considered new point.
!   No ponderation is used (see grille_p)
!
!           (c)
!        ----d-----
!        | . . . .|
!        |        |
!     (b)a . * . .b(a)
!        |        |
!        | . . . .|
!        ----c-----
!           (d)
!
!-------------------------------------------------------------------------------
  IMPLICIT NONE
!-------------------------------------------------------------------------------
! Arguments:
  REAL, INTENT(IN)  :: xdata(:), ydata(:)       !--- INPUT  FIELD X AND Y COORD.
  REAL, INTENT(IN)  :: entree(size(xdata),size(ydata)) !--- INPUT  FIELD
  REAL, INTENT(IN)  :: x(:), y(:)               !--- OUTPUT FIELD X AND Y COORD.
  REAL, INTENT(OUT) :: sortie(size(x),size(y))  !--- OUTPUT FIELD 
!-------------------------------------------------------------------------------
  CALL fine2coarse(xdata,ydata,x,y,sortie,entree)

END SUBROUTINE grille_m
!
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!
SUBROUTINE rugosite(xdata, ydata, entree, x, y, sortie, mask)
!
!-------------------------------------------------------------------------------
! Author:  Z.X. Li (april 1st 1994)
!-------------------------------------------------------------------------------
! Purpose: Remap rugosity length ; constant value (0.001) on oceans.
! Naive method  (see grille_m)
!-------------------------------------------------------------------------------
  IMPLICIT NONE
!-------------------------------------------------------------------------------
! Arguments:
  REAL, INTENT(IN)  :: xdata(:), ydata(:)      !--- INPUT  FIELD X AND Y COORD.
  REAL, INTENT(IN)  :: entree(size(xdata),size(ydata)) !--- INPUT  FIELD
  REAL, INTENT(IN)  :: x(:), y(:)              !--- OUTPUT FIELD X AND Y COORD.
  REAL, INTENT(OUT) :: sortie(size(x),size(y)) !--- OUTPUT FIELD 
  REAL, INTENT(IN)  :: mask  (size(x),size(y)) !--- MASK
!-------------------------------------------------------------------------------
  CALL fine2coarse(xdata,ydata,x,y,sortie,log(entree))
  WHERE(nint(mask)==1)
    sortie(:,:)=exp(sortie(:,:))
  ELSE WHERE
    sortie(:,:)=0.001
  END WHERE

END SUBROUTINE rugosite
!
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!
SUBROUTINE sea_ice(xdata, ydata, glace01, x, y, frac_ice)
!
!-------------------------------------------------------------------------------
! Author:  Z.X. Li (april 1st 1994)
! Purpose: Regrid ice indicator (0 or 1) on a coarser grid to get an ice fract.
! field (between 0 and 1).
! Naive method  (see grille_m)
!-------------------------------------------------------------------------------
  IMPLICIT NONE
!-------------------------------------------------------------------------------
! Arguments:
  REAL, INTENT(IN)  :: xdata(:), ydata(:)      !--- INPUT  FIELD X AND Y COORD.
  REAL, INTENT(IN)  :: glace01(:,:)            !--- INPUT  FIELD
  REAL, INTENT(IN)  :: x(:), y(:)              !--- OUTPUT FIELD X AND Y COORD.
  REAL, INTENT(OUT) :: frac_ice(size(x),size(y)) !--- OUTPUT FIELD 
!-------------------------------------------------------------------------------
  CALL fine2coarse(xdata,ydata,x,y,frac_ice,msk=nint(glace01)==1)

END SUBROUTINE sea_ice
!
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!
SUBROUTINE rugsoro(xrel, yrel, relief, xmod, ymod, rugs)
!
!-------------------------------------------------------------------------------
! Author:  Z.X. Li (april 1st 1994) ; Modifications: august 23rd 1995.
! Purpose: Compute rugosity due to orography by using standard dev in a 1x1 cell
!-------------------------------------------------------------------------------
  IMPLICIT NONE
!-------------------------------------------------------------------------------
! Arguments:
  REAL, INTENT(IN)  :: xrel(:), yrel(:)        !--- INPUT  FIELD X AND Y COORD.
  REAL, INTENT(IN)  :: relief(:,:)             !--- INPUT  FIELD
  REAL, INTENT(IN)  :: xmod(:), ymod(:)        !--- OUTPUT FIELD X AND Y COORD.
  REAL, INTENT(OUT) :: rugs(size(xmod),size(ymod)) !--- OUTPUT FIELD 
!-------------------------------------------------------------------------------
! Local variable:
  INTEGER           :: k, nn
  INTEGER, PARAMETER:: itmp=360, jtmp=180
  REAL  :: out(size(xmod),size(xmod)), amin, amax
  REAL  :: cham1tmp(itmp,jtmp), cham2tmp(itmp,jtmp), xtmp(itmp), ytmp(jtmp)
!-------------------------------------------------------------------------------

!--- TEST INPUT FILE FITS FOR ONE DEGREE RESOLUTION
  xtmp(:)=4.0*atan(1.0)*[(-1.0+REAL(2*k-1)/REAL(itmp),k=1,itmp)]
  ytmp(:)=2.0*atan(1.0)*[(-1.0+REAL(2*k-1)/REAL(jtmp),k=1,jtmp)]
  CALL fine2coarse(xrel,yrel,xtmp,ytmp,cham1tmp,relief,d_o2=cham2tmp)
  cham2tmp(:,:)=cham2tmp(:,:)-cham1tmp(:,:)**2
  nn=count(cham2tmp<=-7.5)
  IF(nn/=0) THEN
    print*,'Problem for rugsoro ; std**2 < -7.5 for several points: ',nn
    CALL abort_gcm("", "", 1)
  END IF
  nn=count(cham2tmp<0)
  IF(nn/=0) print*,'Problem for rugsoro ; std**2 < 0. for several points: ',nn
  WHERE(cham2tmp<0.0) cham2tmp=0.0
  cham2tmp(:,:)=sqrt(cham2tmp(:,:))
  amin=minval(cham2tmp); amax=maxval(cham2tmp)
  print*, 'Ecart-type 1x1:', amin, amax

!--- COMPUTE RUGOSITY AT REQUIRED SCALE
  WHERE(cham2tmp<0.001) cham2tmp=0.001
  CALL fine2coarse(xtmp,ytmp,xmod,ymod,out,REAL(log(cham2tmp)))
  out=exp(out)
  amin=minval(out); amax=maxval(out)
  print*, 'Ecart-type du modele:', amin, amax
  out=out/amax*20.0
  amin=minval(out); amax=maxval(out)
  print*, 'Longueur de rugosite du modele:', amin, amax
  rugs=REAL(out)

END SUBROUTINE rugsoro
!
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!
SUBROUTINE dist_sphe(rf_lon,rf_lat,rlon,rlat,distance)
!
!-------------------------------------------------------------------------------
! Author:  Laurent Li (december 30th 1996).
! Purpose: Compute min. distance (along big circle) between 2 points in radians.
!-------------------------------------------------------------------------------
  IMPLICIT NONE
!-------------------------------------------------------------------------------
! Arguments:
  REAL, INTENT(IN) :: rf_lon, rf_lat  !--- Reference point coordinates (radians)
  REAL, INTENT(IN) :: rlon(:), rlat(:)!--- Points longitudes/latitudes (radians)
  REAL, INTENT(OUT):: distance(size(rlon),size(rlat)) !--- Distance    (radians)
!-------------------------------------------------------------------------------
! Local variables:
  LOGICAL, SAVE :: first=.true.
  REAL,    SAVE :: pi, hpi
  REAL    :: pa, pb, cpa, spa, crlo(size(rlon))
  INTEGER :: i, j
!-------------------------------------------------------------------------------
  IF(first) then; pi=4.0*atan(1.0); hpi=pi/2.0; first=.false.; END IF
  crlo(:)=cos(rf_lon-rlon(:))         !--- COS of points 1 and 2 angle
  pa=hpi-rf_lat                       !--- North Pole - Point 1 distance
  cpa=cos(pa); spa=sin(pa)
  DO j=1,SIZE(rlat)
    pb=hpi-rlat(j)                    !--- North Pole - Point 2 distance
    distance(:,j)=acos(cpa*cos(pb)+spa*sin(pb)*crlo(:))
  END DO

END SUBROUTINE dist_sphe
!
!-------------------------------------------------------------------------------

END MODULE grid_atob_m
!
!*******************************************************************************