Skip to content
Snippets Groups Projects
icedyn_rhg_evp.F90 68.4 KiB
Newer Older
Guillaume Samson's avatar
Guillaume Samson committed
            ztmp(:) = glob_sum_vec( 'icedyn_rhg_evp', zres )
            IF( ztmp(2) /= 0._wp )   zresm = ztmp(1) / ztmp(2)
         ENDIF
      ENDIF

      IF( lwm ) THEN
         ! write variables
         istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) )
         ! close file
         IF( kt == nitend - nn_fsbc + 1 .AND. kiter == kitermax )   istatus = NF90_CLOSE(ncvgid)
      ENDIF

   END SUBROUTINE rhg_cvg


   SUBROUTINE rhg_evp_rst( cdrw, kt )
      !!---------------------------------------------------------------------
      !!                   ***  ROUTINE rhg_evp_rst  ***
      !!
      !! ** Purpose :   Read or write RHG file in restart file
      !!
      !! ** Method  :   use of IOM library
      !!----------------------------------------------------------------------
      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
      !
      INTEGER  ::   iter            ! local integer
      INTEGER  ::   id1, id2, id3   ! local integers
      !!----------------------------------------------------------------------
      !
      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
         !                                   ! ---------------
         IF( ln_rstart ) THEN                   !* Read the restart file
            !
            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
            !
            IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist
               CALL iom_get( numrir, jpdom_auto, 'stress1_i' , stress1_i , cd_type = 'T' )
               CALL iom_get( numrir, jpdom_auto, 'stress2_i' , stress2_i , cd_type = 'T' )
               CALL iom_get( numrir, jpdom_auto, 'stress12_i', stress12_i, cd_type = 'F' )
            ELSE                                     ! start rheology from rest
               IF(lwp) WRITE(numout,*)
               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without rheology, set stresses to 0'
               stress1_i (:,:) = 0._wp
               stress2_i (:,:) = 0._wp
               stress12_i(:,:) = 0._wp
            ENDIF
         ELSE                                   !* Start from rest
            IF(lwp) WRITE(numout,*)
            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set stresses to 0'
            stress1_i (:,:) = 0._wp
            stress2_i (:,:) = 0._wp
            stress12_i(:,:) = 0._wp
         ENDIF
         !
      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
         !                                   ! -------------------
         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
         !
         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i )
         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i )
         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
         !
      ENDIF
      !
   END SUBROUTINE rhg_evp_rst

   SUBROUTINE rhg_upstream( jter, pdt, pu, pv, pt )
      !!---------------------------------------------------------------------
      !!                    ***  ROUTINE rhg_upstream  ***
      !!
      !! **  Purpose :   compute the upstream fluxes and upstream guess of tracer
      !!----------------------------------------------------------------------
      INTEGER                    , INTENT(in   ) ::   jter
      REAL(wp)                   , INTENT(in   ) ::   pdt              ! tracer time-step
      REAL(wp), DIMENSION(:,:  ) , INTENT(in   ) ::   pu, pv           ! 2 ice velocity components
      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   pt               ! tracer fields
      !
      INTEGER  ::   ji, jj, jl    ! dummy loop indices
      REAL(wp) ::   ztra          ! local scalar
      LOGICAL  ::   ll_upsxy = .TRUE.
      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfv_ups, zpt   ! upstream fluxes and tracer guess
      !!----------------------------------------------------------------------
      DO jl = 1, jpl
         IF( .NOT. ll_upsxy ) THEN         !** no alternate directions **!
            !
            DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
               zfu_ups(ji,jj) = MAX(pu(ji,jj)*e2u(ji,jj), 0._wp) * pt(ji,jj,jl) + MIN(pu(ji,jj)*e2u(ji,jj), 0._wp) * pt(ji+1,jj,jl)
               zfv_ups(ji,jj) = MAX(pv(ji,jj)*e1v(ji,jj), 0._wp) * pt(ji,jj,jl) + MIN(pv(ji,jj)*e1v(ji,jj), 0._wp) * pt(ji,jj+1,jl)
            END_2D
            !
         ELSE                              !** alternate directions **!
            !
            IF( MOD(jter,2) == 1 ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
               !
               DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls )       !-- flux in x-direction
                  zfu_ups(ji,jj) = MAX( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji  ,jj,jl) + &
                     &             MIN( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji+1,jj,jl)
               END_2D
               !
               DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls )     !-- first guess of tracer from u-flux
                  ztra       = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) )
                  zpt(ji,jj) =   ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
               END_2D
               !
               DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 )   !-- flux in y-direction
                  zfv_ups(ji,jj) = MAX( pv(ji,jj)*e1v(ji,jj), 0._wp ) * zpt(ji,jj  ) + &
                     &             MIN( pv(ji,jj)*e1v(ji,jj), 0._wp ) * zpt(ji,jj+1)
               END_2D
               !
            ELSE                          !==  even ice time step:  adv_y then adv_x  ==!
               !
               DO_2D( nn_hls, nn_hls, nn_hls, nn_hls-1 )       !-- flux in y-direction
                  zfv_ups(ji,jj) = MAX( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji,jj  ,jl) + &
                     &             MIN( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji,jj+1,jl)
               END_2D
               !
               DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls-1 )     !-- first guess of tracer from v-flux
                  ztra       = - ( zfv_ups(ji,jj) - zfv_ups(ji,jj-1) )
                  zpt(ji,jj) =   ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
               END_2D
               !
               DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 )   !-- flux in x-direction
                  zfu_ups(ji,jj) = MAX( pu(ji,jj)*e2u(ji,jj), 0._wp ) * zpt(ji  ,jj) + &
                     &             MIN( pu(ji,jj)*e2u(ji,jj), 0._wp ) * zpt(ji+1,jj)
               END_2D
               !
            ENDIF
            !
         ENDIF
         !
         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
            ztra         = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) + zfv_ups(ji,jj) - zfv_ups(ji,jj-1) )
            pt(ji,jj,jl) =   ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1)
         END_2D
      END DO
      !
   END SUBROUTINE rhg_upstream

#else
   !!----------------------------------------------------------------------
   !!   Default option         Empty module           NO SI3 sea-ice model
   !!----------------------------------------------------------------------
#endif

   !!==============================================================================
END MODULE icedyn_rhg_evp