source: TOOLS/MOSAIX/src/MOSAIX/interpol.f90 @ 3622

! -*- Mode: f90 -*-
PROGRAM interpol
!!!
!!! Generates interpolation weights from source (src) grid to destination (dst) grid
!!! Most of the configuration is done by the xml file
!!!
   USE xios
   USE mod_wait
   USE mpi
   IMPLICIT NONE
   ! INCLUDE "mpif.h"
   INTEGER, PARAMETER :: rl = SELECTED_REAL_KIND(15,307) !< Default real precision (double)
   INTEGER :: ierr, comm, rank

   TYPE (xios_context) :: ctx_hdl
   INTEGER :: ni_src, nj_src !< Dimensions of the source grid
   INTEGER :: ni_dst, nj_dst !< Dimensions of the source grid
   REAL (kind=rl), ALLOCATABLE :: imask_src (:,:), imask_dst (:,:)
   REAL (kind=rl), ALLOCATABLE :: lon_src(:,:), lat_src(:,:), field_src(:,:)
   LOGICAL, ALLOCATABLE :: lmask_src (:,:), lmask_dst (:,:)
   INTEGER :: nout = 0, jf
   CHARACTER (LEN=20) :: nchar, type_src, type_dst
   LOGICAL :: l_mask_src = .TRUE. , l_mask_dst = .TRUE.
   !
   REAL (kind=rl), PARAMETER :: rpi = ACOS ( -1.0_rl)
   REAL (kind=rl), PARAMETER :: rad = rpi / 180.0_rl
   !
   INTEGER :: narg, ja !< To handle command line arguments
   CHARACTER (LEN=80) :: cmd_arg
   !!
   !< Initialisation
   WRITE (UNIT=nout, FMT="('-- start interpol')" ) 
   CALL MPI_INIT (ierr)

   CALL init_wait
   CALL xios_initialize ("interpol", return_comm=comm)

   CALL MPI_COMM_RANK (comm, rank, ierr)
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- mpi_comm_rank')" ) rank
   !CALL MPI_COMM_SIZE (comm, size, ierr)

   !< Read command line
   narg = iargc ()
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- narg : ', 1I3)" ) rank, narg

   !< Read arguments on command line
   DO ja = 1, narg
      CALL getarg ( ja, cmd_arg )
      SELECT CASE ( TRIM (cmd_arg) )
      CASE ( '--mask_src=yes' ) ; l_mask_src = .TRUE.
      CASE ( '--mask_src=no'  ) ; l_mask_src = .FALSE.
      CASE ( '--mask_dst=yes' ) ; l_mask_dst = .TRUE.
      CASE ( '--mask_dst=no'  ) ; l_mask_dst = .FALSE.
      END SELECT
   END DO

   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- mask_src : ', 1L)" ) rank, l_mask_src
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- mask_dst : ', 1L)" ) rank, l_mask_dst

   !< Context interpol_read : read masks
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- debut interpol_read')" ) rank
   CALL xios_context_initialize  ("interpol_read", comm)
   CALL xios_get_handle          ("interpol_read", ctx_hdl)
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- set current context interpol_read')" ) rank
   CALL xios_set_current_context (ctx_hdl)
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- close context definition interpol_read')" ) rank
   CALL xios_close_context_definition ()

   !< Read characteristics of the source grid
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- lecture domain_src')" ) rank
   CALL xios_get_domain_attr ("domain_src", ni=ni_src, nj=nj_src)
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- lecture domain_src ', 6I9)") rank, ni_src, nj_src
   ALLOCATE ( lon_src (ni_src, nj_src), lat_src (ni_src, nj_src) )
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- lecture lon lat src ', 6I9)") rank, SIZE (lon_src), SHAPE (lon_src), SIZE (lat_src), SHAPE (lat_src)
   !!
   !
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- lecture type src ', 6I9)") rank
   CALL xios_get_domain_attr ("domain_src", TYPE=type_src )
   SELECT CASE ( TRIM (type_src))
   CASE ( "rectilinear" )
      WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- lecture lon lat src rectilinear', 6I9)") rank
      CALL xios_get_domain_attr ("domain_src", lonvalue_1d=lon_src(:,1), latvalue_1d=lat_src(1,:) )
      lon_src (:,:) = SPREAD ( lon_src(:,1), DIM=2, ncopies=nj_src)
      lat_src (:,:) = SPREAD ( lat_src(1,:), DIM=1, ncopies=nj_src)
   CASE default
      IF ( nj_src == 1 ) THEN
         WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- lecture lon lat src 1D ', 6I9)") rank
         CALL xios_get_domain_attr ("domain_src", lonvalue_1d=lon_src(:,1), latvalue_1d=lat_src(:,1) )
      ELSE
         WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- lecture lon lat src 2D ', 6I9)") rank
         CALL xios_get_domain_attr ("domain_src", lonvalue_2d=lon_src(:,:), latvalue_2d=lat_src(:,:) )
      ENDIF
   END SELECT
   !< Read mask on the source grid
   ALLOCATE ( imask_src (ni_src, nj_src), lmask_src (ni_src, nj_src) )
   IF ( l_mask_src ) THEN 
      WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- receive field mask_src')" ) rank
      CALL xios_recv_field ("mask_src", imask_src)
      lmask_src = .FALSE.
      WHERE (imask_src>0) lmask_src = .TRUE.
   ELSE
      imask_src (:,:) = 1 ; lmask_src (:,:) = .TRUE.
   ENDIF
   
   !< Read mask on the destination grid
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- lecture domain_dst')" ) rank
   CALL xios_get_domain_attr ("domain_dst", ni=ni_dst, nj=nj_dst)
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- lecture domain_dst ', 6I7)") rank, ni_dst, nj_dst
   ALLOCATE ( imask_dst (ni_dst, nj_dst), lmask_dst (ni_dst, nj_dst) )
   IF ( l_mask_dst ) THEN
      WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- receive field mask_dst')" ) rank 
      CALL xios_recv_field ("mask_dst", imask_dst)
      lmask_dst = .FALSE.
      WHERE (imask_dst>0) lmask_dst = .TRUE.
   ELSE
      imask_dst (:,:) = 1 ; lmask_dst (:,:) = .TRUE.
   ENDIF

   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- fin interpol_read (1)')") rank
   CALL xios_context_finalize ()
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- fin interpol_read (2)')") rank
   
   !< Context interpol run : generates weights, interpolate mask from source to destination
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- debut interpol_run')" ) rank
   CALL xios_context_initialize  ("interpol_run", comm)
   CALL xios_get_handle          ("interpol_run", ctx_hdl)
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- set context interpol_run')" ) rank
   CALL xios_set_current_context (ctx_hdl)

   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- get attributes domain_src')" ) rank
   IF ( l_mask_src ) CALL xios_set_domain_attr ("domain_src", mask_2d=lmask_src)
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- get attributes domain_src')" ) rank
   IF ( l_mask_dst ) CALL xios_set_domain_attr ("domain_dst", mask_2d=lmask_dst)

   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- close context definition')" ) rank
   CALL xios_close_context_definition ()

   CALL xios_update_calendar (1)
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- send field mask_src')" ) rank
   CALL xios_send_field ("mask_src", imask_src)

   !!< Creates analytic fields and interpolate
   ALLOCATE ( field_src (ni_src, nj_src) )

   DO jf = 1, 6
      WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- working on test case ', 1I2.2)" ) rank, jf
      SELECT CASE ( jf)
      CASE ( 1) ; field_src (:,:) = REAL ( imask_src, kind=rl)
      CASE ( 2) ; field_src (:,:) = SIN ( rad * lat_src(:,:) )
      CASE ( 3) ; field_src (:,:) = COS ( rad * lat_src(:,:) )
      CASE ( 4) ; field_src (:,:) = COS ( rad * lat_src(:,:) ) * COS ( rad * lon_src(:,:) )
      CASE ( 5) ; field_src (:,:) = SIN ( rad * lon_src(:,:) )
      CASE ( 6) ; field_src (:,:) = SIN ( rad * lon_src(:,:) + rpi/2.0_rl )
      END SELECT
      WRITE (UNIT=nchar, FMT='("field", 1I2.2, "_src")' ) jf
      CALL xios_send_field ( TRIM(nchar), field_src )
   END DO
   
   !!<
   
   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- fin interpol_run')" ) rank
   CALL xios_context_finalize ()

   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- mpi_comm_free')" ) rank
   CALL MPI_COMM_FREE (comm, ierr)

   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- xios finalize')" ) rank
   CALL xios_finalize ()

   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- mpi finalize')" ) rank
   CALL MPI_FINALIZE (ierr)

   WRITE (UNIT=nout, FMT="('-- ', 1I4.4, ' -- fini')" ) rank

END PROGRAM interpol