write_data Subroutine

public subroutine write_data(this, rc)

Type Bound

RegridSupport

Arguments

Type IntentOptional Attributes Name
class(RegridSupport), intent(inout), target :: this
integer, intent(out), optional :: rc

Calls

proc~~write_data~~CallsGraph proc~write_data RegridSupport%write_data ESMF_FieldRegridStore ESMF_FieldRegridStore proc~write_data->ESMF_FieldRegridStore ESMF_GridGetCoord ESMF_GridGetCoord proc~write_data->ESMF_GridGetCoord ESMF_VMBarrier ESMF_VMBarrier proc~write_data->ESMF_VMBarrier ESMF_VMGetGlobal ESMF_VMGetGlobal proc~write_data->ESMF_VMGetGlobal at at proc~write_data->at esmf_fieldcreate esmf_fieldcreate proc~write_data->esmf_fieldcreate get_index get_index proc~write_data->get_index lev lev proc~write_data->lev none~first~177 StringVariableMapIterator%first proc~write_data->none~first~177 none~ftn_begin~25 StringVariableMap%ftn_begin proc~write_data->none~ftn_begin~25 none~ftn_end~25 StringVariableMap%ftn_end proc~write_data->none~ftn_end~25 none~get_attribute Variable%get_attribute proc~write_data->none~get_attribute none~get_dimensions Variable%get_dimensions proc~write_data->none~get_dimensions none~get_values UnlimitedEntity%get_values proc~write_data->none~get_values none~get_variable FileMetadata%get_variable proc~write_data->none~get_variable none~get_variables FileMetadata%get_variables proc~write_data->none~get_variables none~get_var~21 NetCDF4_FileFormatter%get_var proc~write_data->none~get_var~21 none~next~82 StringVariableMapIterator%next proc~write_data->none~next~82 none~put_var~21 NetCDF4_FileFormatter%put_var proc~write_data->none~put_var~21 proc~cosd cosd proc~write_data->proc~cosd proc~mapl_verify MAPL_Verify proc~write_data->proc~mapl_verify proc~regrid regrid proc~write_data->proc~regrid proc~sind sind proc~write_data->proc~sind

Called by

proc~~write_data~~CalledByGraph proc~write_data RegridSupport%write_data program~main~2 main program~main~2->proc~write_data

Source Code

   subroutine write_data(this, rc)
      class (RegridSupport), target, intent(inout) :: this
      integer, optional, intent(out) :: rc


      type (StringVariableMapIterator) :: var_iter
      type (StringVariableMap), pointer :: variables

      real(kind=REAL32), target, allocatable :: cs_scalar_patch(:,:)
      real(kind=REAL32), target, allocatable :: cs_vector_patch(:,:,:)
      real(kind=REAL32), target, allocatable :: cs_uvw(:,:,:)
      real(kind=REAL32), target, allocatable :: ll_uvw(:,:,:)
      real(kind=REAL32), target, allocatable :: ll_scalar_patch(:,:)
      real(kind=REAL32), target, allocatable :: ll_vector_patch(:,:,:)

      character(len=:), pointer :: var_name
      type (ESMF_Field) :: cs_scalar_field
      type (ESMF_Field) :: ll_scalar_field
      type (ESMF_Field) :: cs_uvw_field(3)
      type (ESMF_Field) :: ll_uvw_field(3)
      integer :: level
      integer :: num_levels
      integer :: status
      integer :: time
      type (StringVector), pointer :: dims
      type (Variable), pointer :: var

      integer, allocatable :: cs_start(:), cs_count(:)
      integer, allocatable :: ll_start(:), ll_count(:)

      real(kind=REAL64), pointer :: lat(:,:)
      real(kind=REAL64), pointer :: lon(:,:)

      type (Attribute), pointer :: missing_attr
      class(*), pointer :: missing_ptr(:)
      real(kind=REAL32), pointer :: missing

      integer :: d
      logical :: is_scalar
      logical :: is_east_vector_component
      integer :: idx
      character(len=:), allocatable :: north_component
      integer(kind=INT64) :: c0, c1,crate

      associate (cs_fmtr => this%formatter_cubed_sphere, ll_fmtr => this%formatter_lat_lon)
      call cs_fmtr%open(this%in_file, mode=pFIO_READ, rc=status)
      if (status /= pFIO_SUCCESS) then
         if (local_pet == 0) then
            print*, 'Unable to open input file: ',this%in_file
          end if
      end if
      call ll_fmtr%open(this%out_file, mode=pFIO_WRITE, rc=status)
      if (status /= pFIO_SUCCESS) then
         if (local_pet == 0) then
            print*, 'Unable to open new output file: ',this%out_file
            print*, 'Possibly it already exists?'
         end if
      end if

      call ll_fmtr%put_var('lat', this%latitudes)
      call ll_fmtr%put_var('lon', this%longitudes)

      block
        real(kind=REAL64) :: lev(this%LM)
        if (associated(this%cfio_cubed_sphere%get_variable('lev'))) then
           call cs_fmtr%get_var('lev', lev)
           call ll_fmtr%put_var('lev', lev)
        end if
      end block
      block
        integer :: time(this%NT)
        call cs_fmtr%get_var('time', time)
        call ll_fmtr%put_var('time', time)
      end block

      allocate(ll_scalar_patch(this%i_1:this%i_n, this%j_1:this%j_n))
      allocate(ll_vector_patch(this%i_1:this%i_n, this%j_1:this%j_n,2))
      allocate(ll_uvw(this%i_1:this%i_n, this%j_1:this%j_n,3))


      allocate(cs_scalar_patch(this%nx_loc,this%ny_loc))
      allocate(cs_vector_patch(this%nx_loc,this%ny_loc,2))

      allocate(cs_uvw(this%nx_loc,this%ny_loc,3))

      ! Create fields
      ll_scalar_field = ESMF_FieldCreate(this%grid_lat_lon, ll_scalar_patch(:,:), ESMF_INDEX_DELOCAL, &
           & datacopyflag=ESMF_DATACOPY_REFERENCE, rc=status)
      _VERIFY(status)
      cs_scalar_field = ESMF_FieldCreate(this%grid_cubed_sphere, cs_scalar_patch(:,:), ESMF_INDEX_DELOCAL, &
           & datacopyflag=ESMF_DATACOPY_REFERENCE, rc=status)
      _VERIFY(status)
      do d = 1, 3
         ll_uvw_field(d) = ESMF_FieldCreate(this%grid_lat_lon, ll_uvw(:,:,d), ESMF_INDEX_DELOCAL, &
              & datacopyflag=ESMF_DATACOPY_REFERENCE, rc=status)
         _VERIFY(status)
         cs_uvw_field(d) = ESMF_FieldCreate(this%grid_cubed_sphere, cs_uvw(:,:,d), ESMF_INDEX_DELOCAL, &
              & datacopyflag=ESMF_DATACOPY_REFERENCE, rc=status)
         _VERIFY(status)
      end do

      ! Create regrid
      srcTerm = 1
      call system_clock(c0,crate)
      call ESMF_FieldRegridStore(cs_scalar_field, ll_scalar_field, &
           & regridmethod=ESMF_REGRIDMETHOD_BILINEAR, lineType=ESMF_LINETYPE_GREAT_CIRCLE, &
           & srcTermProcessing=srcTerm, &
           & routehandle=default_route_handle, rc=status)
      _VERIFY(status)
      call system_clock(c1,crate)
      if (local_pet == 0) then
         print*,'regrid store: ', real(c1-c0)/crate
      end if

      block
        type (ESMF_VM) :: global
        call ESMF_VMGetGlobal(global, rc=status)
        call ESMF_VMBarrier(global, rc=status)
      end block
      variables => this%cfio_cubed_sphere%get_variables()
      var_iter = variables%ftn_begin()
      do while (var_iter /= variables%ftn_end())
         call var_iter%next()

         var_name => var_iter%first()

         select case (var_name)
         case ('nf', 'ncontact', 'cubed_sphere', &
              & 'Xdim', 'Ydim', 'lons', 'lats', &
              & 'contacts', 'orientation', 'anchor', &
              & 'lev','time')
            ! skip
         case default

            if ( local_pet == 0 .and. this%debug) then
               print*, 'var = ', var_name
            end if

            var => var_iter%second()
            missing_attr => var%get_attribute('missing_value')
            missing_ptr => missing_attr%get_values()

            select type (missing_ptr)
            type is (real(kind=REAL64))
               allocate(missing) ! memory leak!
               missing = missing_ptr(1)
            type is (real(kind=REAL32))
               missing => missing_ptr(1)
            class default
               ! no missing value?
            end select
            num_levels = 1
            dims => var%get_dimensions()
            if (dims%get_index('lev') /= 0) then
               num_levels = this%LM
            else
               num_levels = 1
            end if

            is_scalar = .false.
            is_east_vector_component = .false.

            do idx = 1, this%scalar_variables%size()
               if (this%scalar_variables%at(idx) == var_name) then
                  is_scalar = .true.
                  exit
               end if
            end do

            if (.not. is_scalar) then
               do idx = 1, this%vector_variables(1)%size()
                  if (this%vector_variables(1)%at(idx) == var_name) then
                     is_east_vector_component = .true.
                     north_component = this%vector_variables(2)%at(idx)
                     exit
                  end if
               end do
            end if

            if (.not. (is_scalar .or. is_east_vector_component)) then
               cycle
            end if


            do time = 1, this%nt
               do level = 1, num_levels

                  if (num_levels == 1) then
                     cs_start = [this%x_1,this%y_1,this%my_tile,time]
                     ll_start = [this%i_1,this%j_1,time]
                     cs_count = [this%nx_loc,this%ny_loc,1,1]
                     ll_count = [this%i_n-this%i_1+1, this%j_n-this%j_1+1,1]
                  else
                     cs_start = [this%x_1,this%y_1,this%my_tile,level,time]
                     ll_start = [this%i_1,this%j_1,level,time]
                     cs_count = [this%nx_loc,this%ny_loc,1,1,1]
                     ll_count = [this%i_n-this%i_1+1, this%j_n-this%j_1+1,1,1]
                  end if

                  if (is_east_vector_component) then ! vector

                     call cs_fmtr%get_var(var_name, cs_vector_patch(:,:,1), start=cs_start,count=cs_count, rc=status)
                     _VERIFY(status)
                     call cs_fmtr%get_var(north_component, cs_vector_patch(:,:,2), start=cs_start,count=cs_count, rc=status)
                     _VERIFY(status)

                     call ESMF_GridGetCoord(this%grid_cubed_sphere, coordDim=1, localDE=0, &
                          & staggerLoc=ESMF_STAGGERLOC_CENTER, farrayPtr=lon,rc=status)
                     call ESMF_GridGetCoord(this%grid_cubed_sphere, coordDim=2, localDE=0, &
                          & staggerLoc=ESMF_STAGGERLOC_CENTER, farrayPtr=lat,rc=status)

                     associate (u => cs_vector_patch(:,:,1), v => cs_vector_patch(:,:,2))
                       where (u == missing)
                          cs_uvw(:,:,1) = missing
                          cs_uvw(:,:,2) = missing
                          cs_uvw(:,:,3) = missing
                       elsewhere
                          cs_uvw(:,:,1) = - u * sind(lon) - v * sind(lat)*cosd(lon)
                          cs_uvw(:,:,2) = + u * cosd(lon) - v * sind(lat)*sind(lon)
                          cs_uvw(:,:,3) = + v * cosd(lat)
                      end where
                     end associate


                     do d = 1, 3
                        if (associated(missing)) then
                           call regrid(srcField=cs_uvw_field(d), dstField=ll_uvw_field(d), missing=missing, rc=status)
                           _VERIFY(status)
                        else
                           call regrid(srcField=cs_uvw_field(d), dstField=ll_uvw_field(d), rc=status)
                           _VERIFY(status)
                        end if
                     end do

                     call ESMF_GridGetCoord(this%grid_lat_lon, coordDim=1, localDE=0, &
                          & staggerLoc=ESMF_STAGGERLOC_CENTER, farrayPtr=lon,rc=status)
                     call ESMF_GridGetCoord(this%grid_lat_lon, coordDim=2, localDE=0, &
                          & staggerLoc=ESMF_STAGGERLOC_CENTER, farrayPtr=lat,rc=status)

                     associate (u => ll_vector_patch(:,:,1), v => ll_vector_patch(:,:,2))
                       where (ll_uvw(:,:,1) == missing)
                          u = missing
                          v = missing
                       elsewhere
                          u = -ll_uvw(:,:,1) * sind(lon) + ll_uvw(:,:,2) * cosd(lon)
                          v = -ll_uvw(:,:,1) * sind(lat)*cosd(lon) - ll_uvw(:,:,2) * sind(lat)*sind(lon) +  ll_uvw(:,:,3)*cosd(lat)
                       end where
                     end associate

                     call ll_fmtr%put_var(var_name, ll_vector_patch(:,:,1), start=ll_start, count=ll_count, rc=status)
                     _VERIFY(status)
                     call ll_fmtr%put_var(north_component, ll_vector_patch(:,:,2), start=ll_start, count=ll_count, rc=status)
                     _VERIFY(status)

                  else ! scalar
                     call cs_fmtr%get_var(var_name, cs_scalar_patch, start=cs_start,count=cs_count, rc=status)
                     _VERIFY(status)
                     if (associated(missing)) then
                        call regrid(srcField=cs_scalar_field, dstField=ll_scalar_field, missing=missing, rc=status)
                        _VERIFY(status)
                     else
                        call regrid(srcField=cs_scalar_field, dstField=ll_scalar_field, rc=status)
                        _VERIFY(status)
                     end if
                     call ll_fmtr%put_var(var_name, ll_scalar_patch, start=ll_start, count=ll_count, rc=status)
                     _VERIFY(status)

                  end if

               end do
            end do
         end select
       end do


      call ll_fmtr%close()
      call cs_fmtr%close()

      end associate

   end subroutine write_data