forcing-perturbation/interpolation/rotateUVorca_sources/rotateUVorca.f90

program rotateUVorca
!==============================================================================
! This program rotates U and V components from the geographical directions 
! toward the spherical grid directions based on NEMO3.2 routines 
!
! Written on 2012/02/21
! Author : Virginie Guemas
!==============================================================================
!  
     use par_kind
     use netcdf
     use geo2ocean
     use dom_oce
     use handerr

     implicit none

     include 'netcdf.inc'
  
     character (80) :: &
     & Ufilein,  &  ! filename containing the Eastward component
     & Uvarin,   &  ! name of the Eastward component
     & Vfilein,  &  ! filename containing the Northward component
     & Vvarin,   &  ! name of the Northward component
     & meshmask, &  ! name of the meshmask 
     & Ufileout, &   ! U output file
     & Vfileout     ! V output file

     integer :: nc_fileU_id, nc_fileV_id, nc_filemask_id, nc_varU_id, &
     & nc_varV_id, nc_time_id, nc_var_type, nc_outfile_id, ncstat,    &
     & nc_glamt_id, nc_glamu_id, nc_glamv_id, nc_glamf_id,            &
     & nc_gphit_id, nc_gphiu_id, nc_gphiv_id, nc_gphif_id

     integer, dimension(:), allocatable :: nc_dims_ids
 
     character (80) :: timename

     integer :: ndims, ntime, jtime

     real (kind=wp), dimension(:), allocatable ::           &
     &    time

     real (kind=wp), dimension(:,:,:), allocatable ::       &
     &    Ufield, Vfield, Ufield2, Vfield2
!
!==============================================================================
! 
     namelist /nam_rotUV/ Ufilein, Uvarin, Vfilein, Vvarin, &
     &                    meshmask, Ufileout, Vfileout
!
!==============================================================================
!
!                             Read namelist
!
!============================================================================== 
!
      open(80, file='namelist_rotateUVorca', status='old', form='formatted')
      read(80, nml=nam_rotUV)
      write(*,nml=nam_rotUV)
!
!============================================================================== 
!
!                       Read input (U,V) components
!
!==============================================================================
!
     ncstat = nf_open(Ufilein, NF90_NOWRITE, nc_fileU_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Open file Ufilein   ")
     ncstat = nf_inq_varid(nc_fileU_id, Uvarin, nc_varU_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Inq ID varU         ")
     ncstat = nf_inq_varndims(nc_fileU_id, nc_varU_id, ndims)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Inq dims varU       ")
     if ( ndims < 2 .or. ndims > 3) then
        stop "Input files should have (lon, lat) ot (lon, lat, time) &
     &       dimensions"
     endif

     if ( ndims == 3) then
       allocate(nc_dims_ids(ndims))
       ncstat = nf_inq_vardimid(nc_fileU_id, nc_varU_id, nc_dims_ids)
       if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Inq ID dims varU    ")
       ncstat = nf90_inquire_dimension(nc_fileU_id, nc_dims_ids(3), & 
     & timename, ntime)
       if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Inq. dims time      ")
       allocate(time(ntime))
       time=0.
       !ncstat = nf90_get_var(nc_fileU_id, nc_time_id, time)
       !if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Get var time        ")
     else 
       ntime=1
       allocate(time(1))
       time=1.
     endif

     print*, jpi, jpj, ntime
     allocate( Ufield (jpi, jpj, ntime))
     ncstat = nf90_get_var(nc_fileU_id, nc_varU_id, Ufield)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Get var varU        ")
     ncstat = nf_close(nc_fileU_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Close fileU         ")


     ncstat = nf_open(Vfilein, NF90_NOWRITE, nc_fileV_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Open fileV          ")
     ncstat = nf_inq_varid(nc_fileV_id, Vvarin, nc_varV_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Inq. ID varV        ")
     ncstat = nf_inq_varndims(nc_fileV_id, nc_varV_id, ndims) 
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Inq. dims varV      ")
     if ( ndims /= size(nc_dims_ids) ) then
        stop "Input files should have the same dimensions"
      endif

     allocate( Vfield (jpi, jpj, ntime))
     ncstat = nf90_get_var(nc_fileV_id, nc_varV_id, Vfield)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Get var varV        ")
     ncstat = nf_close(nc_fileV_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Close fileV         ")
!==============================================================================
!
!                           Read meshmask
!
!==============================================================================
!
     ncstat = nf_open(meshmask, NF90_NOWRITE, nc_filemask_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Open filemask       ")

     ncstat = nf_inq_varid(nc_filemask_id, 'glamt', nc_glamt_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Inq ID glamt        ")
     ncstat = nf90_get_var(nc_filemask_id, nc_glamt_id, glamt)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Get var glamt       ")

     ncstat = nf_inq_varid(nc_filemask_id, 'glamf', nc_glamf_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Inq ID glamf        ")
     ncstat = nf90_get_var(nc_filemask_id, nc_glamf_id, glamf)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Get var glamf       ")
  
     ncstat = nf_inq_varid(nc_filemask_id, 'glamu', nc_glamu_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Inq ID glamu        ")
     ncstat = nf90_get_var(nc_filemask_id, nc_glamu_id, glamu)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Get var glamu       ")

     ncstat = nf_inq_varid(nc_filemask_id, 'glamv', nc_glamv_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Inq ID glamv        ")
     ncstat = nf90_get_var(nc_filemask_id, nc_glamv_id, glamv)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Get var glamv       ")

     ncstat = nf_inq_varid(nc_filemask_id, 'gphit', nc_gphit_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Inq ID gphit        ")
     ncstat = nf90_get_var(nc_filemask_id, nc_gphit_id, gphit)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Get var gphit       ")

     ncstat = nf_inq_varid(nc_filemask_id, 'gphif', nc_gphif_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Inq ID gphif        ")
     ncstat = nf90_get_var(nc_filemask_id, nc_gphif_id, gphif)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Get var gphif       ")

     ncstat = nf_inq_varid(nc_filemask_id, 'gphiu', nc_gphiu_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Inq ID gphiu        ")
     ncstat = nf90_get_var(nc_filemask_id, nc_gphiu_id, gphiu)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Get var gphiu       ")

     ncstat = nf_inq_varid(nc_filemask_id, 'gphiv', nc_gphiv_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Inq ID gphiv        ")
     ncstat = nf90_get_var(nc_filemask_id, nc_gphiv_id, gphiv)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Get var gphiv       ")

     ncstat = nf_close(nc_filemask_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat, "Close filemask      ")
!
!==============================================================================
! 
!                          Perform rotation
!
!==============================================================================
!
      allocate( Ufield2 (jpi, jpj, ntime))
      allocate( Vfield2 (jpi, jpj, ntime))
      do jtime = 1,ntime
        call rot_rep(Ufield(:,:,jtime),Vfield(:,:,jtime),'T','en->i',Ufield2(:,:,jtime))
        call rot_rep(Ufield(:,:,jtime),Vfield(:,:,jtime),'T','en->j',Vfield2(:,:,jtime))
      end do
      deallocate(Ufield,Vfield)
!
!============================================================================== 
!
!                        Create output netcdf
!
!============================================================================== 
!
     ncstat = nf_create (Ufileout, NF_CLOBBER, nc_outfile_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat)
     ncstat = nf_def_dim (nc_outfile_id, 'x', jpi, nc_dims_ids(1))
     if (ncstat .ne. nf_noerr) call handle_err(ncstat)
     ncstat = nf_def_dim (nc_outfile_id, 'y', jpj, nc_dims_ids(2))
     if (ncstat .ne. nf_noerr) call handle_err(ncstat)
     if ( ndims > 2) then 
       ncstat = nf_def_dim (nc_outfile_id, 'time', ntime, nc_dims_ids(3))
       if (ncstat .ne. nf_noerr) call handle_err(ncstat)
       ncstat = nf_def_var (nc_outfile_id, 'time', NF_DOUBLE, 1, &
    &           nc_dims_ids(3), nc_time_id)
       if (ncstat .ne. nf_noerr) call handle_err(ncstat)
     endif
     ncstat = nf_def_var (nc_outfile_id, Uvarin, NF_DOUBLE, 3,          &
    &                    nc_dims_ids, nc_varU_id )
     if (ncstat .ne. nf_noerr) call handle_err(ncstat)
     ncstat = nf_enddef(nc_outfile_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat)

     ncstat = nf_put_var_double(nc_outfile_id, nc_varU_id, Ufield2)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat)
     if ( ndims > 2) then
       ncstat = nf_put_var_double(nc_outfile_id, nc_time_id, time)
       if (ncstat .ne. nf_noerr) call handle_err(ncstat)
     endif
     ncstat = nf_close(nc_outfile_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat)
!
!============================================================================== 
!
     ncstat = nf_create (Vfileout, NF_CLOBBER, nc_outfile_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat)
     ncstat = nf_def_dim (nc_outfile_id, 'x', jpi, nc_dims_ids(1))
     if (ncstat .ne. nf_noerr) call handle_err(ncstat)
     ncstat = nf_def_dim (nc_outfile_id, 'y', jpj, nc_dims_ids(2))
     if (ncstat .ne. nf_noerr) call handle_err(ncstat)
     if ( ndims > 2) then 
       ncstat = nf_def_dim (nc_outfile_id, 'time', ntime, nc_dims_ids(3))
       if (ncstat .ne. nf_noerr) call handle_err(ncstat)
       ncstat = nf_def_var (nc_outfile_id, 'time', NF_DOUBLE, 1, &
    &           nc_dims_ids(3), nc_time_id)
       if (ncstat .ne. nf_noerr) call handle_err(ncstat)
     endif
     ncstat = nf_def_var (nc_outfile_id, Vvarin, NF_DOUBLE, 3,          &
    &                    nc_dims_ids, nc_varV_id )
     if (ncstat .ne. nf_noerr) call handle_err(ncstat)
     ncstat = nf_enddef(nc_outfile_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat)

     ncstat = nf_put_var_double(nc_outfile_id, nc_varV_id, Vfield2)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat)
     if ( ndims > 2) then
       ncstat = nf_put_var_double(nc_outfile_id, nc_time_id, time)
       if (ncstat .ne. nf_noerr) call handle_err(ncstat)
     endif
     ncstat = nf_close(nc_outfile_id)
     if (ncstat .ne. nf_noerr) call handle_err(ncstat)
     deallocate(Ufield2,Vfield2,time,nc_dims_ids)

end program rotateUVorca