Skip to content
Snippets Groups Projects
icedyn_adv_pra.F90 72.2 KiB
Newer Older
Guillaume Samson's avatar
Guillaume Samson committed
               CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   syc0 (:,:,jk,:) = z3d(:,:,:)
               znam = 'sxxc0'//'_l'//zchar1
               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxc0(:,:,jk,:) = z3d(:,:,:)
               znam = 'syyc0'//'_l'//zchar1
               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syyc0(:,:,jk,:) = z3d(:,:,:)
               znam = 'sxyc0'//'_l'//zchar1
               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxyc0(:,:,jk,:) = z3d(:,:,:)
            END DO
            !                                                        ! ice layers heat content
            DO jk = 1, nlay_i
               WRITE(zchar1,'(I2.2)') jk
               znam = 'sxe'//'_l'//zchar1
               CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sxe (:,:,jk,:) = z3d(:,:,:)
               znam = 'sye'//'_l'//zchar1
               CALL iom_get( numrir, jpdom_auto, znam , z3d, psgn = -1._wp )   ;   sye (:,:,jk,:) = z3d(:,:,:)
               znam = 'sxxe'//'_l'//zchar1
               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxxe(:,:,jk,:) = z3d(:,:,:)
               znam = 'syye'//'_l'//zchar1
               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   syye(:,:,jk,:) = z3d(:,:,:)
               znam = 'sxye'//'_l'//zchar1
               CALL iom_get( numrir, jpdom_auto, znam , z3d )   ;   sxye(:,:,jk,:) = z3d(:,:,:)
            END DO
            !
            IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN                                    ! melt pond fraction
               IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN
                  CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap , psgn = -1._wp )
                  CALL iom_get( numrir, jpdom_auto, 'syap' , syap , psgn = -1._wp )
                  CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap )
                  CALL iom_get( numrir, jpdom_auto, 'syyap', syyap )
                  CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap )
                  !                                                     ! melt pond volume
                  CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp , psgn = -1._wp )
                  CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp , psgn = -1._wp )
                  CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp )
                  CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp )
                  CALL iom_get( numrir, jpdom_auto, 'sxyvp', sxyvp )
               ELSE
                  sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp   ! melt pond fraction
                  sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp   ! melt pond volume
               ENDIF
                  !
               IF ( ln_pnd_lids ) THEN                               ! melt pond lid volume
                  IF( iom_varid( numrir, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN
                     CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl , psgn = -1._wp )
                     CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl , psgn = -1._wp )
                     CALL iom_get( numrir, jpdom_auto, 'sxxvl', sxxvl )
                     CALL iom_get( numrir, jpdom_auto, 'syyvl', syyvl )
                     CALL iom_get( numrir, jpdom_auto, 'sxyvl', sxyvl )
                  ELSE
                     sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp   ! melt pond lid volume
                  ENDIF
               ENDIF
            ENDIF
            !
         ELSE                                   !**  start rheology from rest  **!
            !
            IF(lwp) WRITE(numout,*) '   ==>>   start from rest OR previous run without Prather, set moments to 0'
            !
            sxice = 0._wp   ;   syice = 0._wp   ;   sxxice = 0._wp   ;   syyice = 0._wp   ;   sxyice = 0._wp      ! ice thickness
            sxsn  = 0._wp   ;   sysn  = 0._wp   ;   sxxsn  = 0._wp   ;   syysn  = 0._wp   ;   sxysn  = 0._wp      ! snow thickness
            sxa   = 0._wp   ;   sya   = 0._wp   ;   sxxa   = 0._wp   ;   syya   = 0._wp   ;   sxya   = 0._wp      ! ice concentration
            sxsal = 0._wp   ;   sysal = 0._wp   ;   sxxsal = 0._wp   ;   syysal = 0._wp   ;   sxysal = 0._wp      ! ice salinity
            sxage = 0._wp   ;   syage = 0._wp   ;   sxxage = 0._wp   ;   syyage = 0._wp   ;   sxyage = 0._wp      ! ice age
            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content
            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content
            IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN
               sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp       ! melt pond fraction
               sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp       ! melt pond volume
               IF ( ln_pnd_lids ) THEN
                  sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp       ! melt pond lid volume
               ENDIF
            ENDIF
         ENDIF
         !
         !                                   !=====================================!
      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   !==  write in the ice restart file  ==!
         !                                   !=====================================!
         IF(lwp) WRITE(numout,*) '----  ice-adv-rst  ----'
         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
         !
         !
         ! In case Prather scheme is used for advection, write second order moments
         ! ------------------------------------------------------------------------
         !
         !                                                           ! ice thickness
         CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice  )
         CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice  )
         CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice )
         CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice )
         CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice )
         !                                                           ! snow thickness
         CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  , sxsn   )
         CALL iom_rstput( iter, nitrst, numriw, 'sysn'  , sysn   )
         CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn  )
         CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn  )
         CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn  )
         !                                                           ! ice concentration
         CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa    )
         CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya    )
         CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  , sxxa   )
         CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya   )
         CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya   )
         !                                                           ! ice salinity
         CALL iom_rstput( iter, nitrst, numriw, 'sxsal' , sxsal  )
         CALL iom_rstput( iter, nitrst, numriw, 'sysal' , sysal  )
         CALL iom_rstput( iter, nitrst, numriw, 'sxxsal', sxxsal )
         CALL iom_rstput( iter, nitrst, numriw, 'syysal', syysal )
         CALL iom_rstput( iter, nitrst, numriw, 'sxysal', sxysal )
         !                                                           ! ice age
         CALL iom_rstput( iter, nitrst, numriw, 'sxage' , sxage  )
         CALL iom_rstput( iter, nitrst, numriw, 'syage' , syage  )
         CALL iom_rstput( iter, nitrst, numriw, 'sxxage', sxxage )
         CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage )
         CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage )
         !                                                           ! snow layers heat content
         DO jk = 1, nlay_s
            WRITE(zchar1,'(I2.2)') jk
            znam = 'sxc0'//'_l'//zchar1  ;   z3d(:,:,:) = sxc0 (:,:,jk,:)
            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
            znam = 'syc0'//'_l'//zchar1  ;   z3d(:,:,:) = syc0 (:,:,jk,:)
            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
            znam = 'sxxc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxxc0(:,:,jk,:)
            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
            znam = 'syyc0'//'_l'//zchar1 ;   z3d(:,:,:) = syyc0(:,:,jk,:)
            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
            znam = 'sxyc0'//'_l'//zchar1 ;   z3d(:,:,:) = sxyc0(:,:,jk,:)
            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
         END DO
         !                                                           ! ice layers heat content
         DO jk = 1, nlay_i
            WRITE(zchar1,'(I2.2)') jk
            znam = 'sxe'//'_l'//zchar1   ;   z3d(:,:,:) = sxe (:,:,jk,:)
            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
            znam = 'sye'//'_l'//zchar1   ;   z3d(:,:,:) = sye (:,:,jk,:)
            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
            znam = 'sxxe'//'_l'//zchar1  ;   z3d(:,:,:) = sxxe(:,:,jk,:)
            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
            znam = 'syye'//'_l'//zchar1  ;   z3d(:,:,:) = syye(:,:,jk,:)
            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
            znam = 'sxye'//'_l'//zchar1  ;   z3d(:,:,:) = sxye(:,:,jk,:)
            CALL iom_rstput( iter, nitrst, numriw, znam , z3d )
         END DO
         !
         IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN                                       ! melt pond fraction
            CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  )
            CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  )
            CALL iom_rstput( iter, nitrst, numriw, 'sxxap', sxxap )
            CALL iom_rstput( iter, nitrst, numriw, 'syyap', syyap )
            CALL iom_rstput( iter, nitrst, numriw, 'sxyap', sxyap )
            !                                                        ! melt pond volume
            CALL iom_rstput( iter, nitrst, numriw, 'sxvp' , sxvp  )
            CALL iom_rstput( iter, nitrst, numriw, 'syvp' , syvp  )
            CALL iom_rstput( iter, nitrst, numriw, 'sxxvp', sxxvp )
            CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp )
            CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp )
            !
            IF ( ln_pnd_lids ) THEN                                  ! melt pond lid volume
               CALL iom_rstput( iter, nitrst, numriw, 'sxvl' , sxvl  )
               CALL iom_rstput( iter, nitrst, numriw, 'syvl' , syvl  )
               CALL iom_rstput( iter, nitrst, numriw, 'sxxvl', sxxvl )
               CALL iom_rstput( iter, nitrst, numriw, 'syyvl', syyvl )
               CALL iom_rstput( iter, nitrst, numriw, 'sxyvl', sxyvl )
            ENDIF
         ENDIF
         !
      ENDIF
      !
   END SUBROUTINE adv_pra_rst

   SUBROUTINE icemax3D( pice , pmax )
      !!---------------------------------------------------------------------
      !!                   ***  ROUTINE icemax3D ***
      !! ** Purpose :  compute the max of the 9 points around
      !!----------------------------------------------------------------------
      REAL(wp), DIMENSION(:,:,:), INTENT(in ) ::   pice   ! input
      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   pmax   ! output
      !
      REAL(wp), DIMENSION(Nis0:Nie0) ::   zmax1, zmax2
      REAL(wp)                       ::   zmax3
      INTEGER  ::   ji, jj, jl   ! dummy loop indices
      !!----------------------------------------------------------------------
      ! basic version: get the max of epsi20 + 9 neighbours
!!$      DO jl = 1, jpl
!!$         DO_2D( 0, 0, 0, 0 )
!!$            pmax(ji,jj,jl) = MAX( epsi20, pice(ji-1,jj-1,jl), pice(ji,jj-1,jl), pice(ji+1,jj-1,jl),   &
!!$               &                          pice(ji-1,jj  ,jl), pice(ji,jj  ,jl), pice(ji+1,jj  ,jl),   &
!!$               &                          pice(ji-1,jj+1,jl), pice(ji,jj+1,jl), pice(ji+1,jj+1,jl) )
!!$         END_2D
!!$      END DO
      ! optimized version : does a little bit more than 2 max of epsi20 + 3 neighbours
      DO jl = 1, jpl
         DO ji = Nis0, Nie0
            zmax1(ji) = MAX( epsi20, pice(ji,Njs0-1,jl), pice(ji-1,Njs0-1,jl), pice(ji+1,Njs0-1,jl) )
            zmax2(ji) = MAX( epsi20, pice(ji,Njs0  ,jl), pice(ji-1,Njs0  ,jl), pice(ji+1,Njs0  ,jl) )
         END DO
         DO_2D( 0, 0, 0, 0 )
            zmax3 = MAX( epsi20, pice(ji,jj+1,jl), pice(ji-1,jj+1,jl), pice(ji+1,jj+1,jl) )
            pmax(ji,jj,jl) = MAX( epsi20, zmax1(ji), zmax2(ji), zmax3 )
            zmax1(ji) = zmax2(ji)
            zmax2(ji) = zmax3
         END_2D
      END DO
   END SUBROUTINE icemax3D

   SUBROUTINE icemax4D( pice , pmax )
      !!---------------------------------------------------------------------
      !!                   ***  ROUTINE icemax4D ***
      !! ** Purpose :  compute the max of the 9 points around
      !!----------------------------------------------------------------------
      REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) ::   pice   ! input
      REAL(wp), DIMENSION(:,:,:,:), INTENT(out) ::   pmax   ! output
      !
      REAL(wp), DIMENSION(Nis0:Nie0) ::   zmax1, zmax2
      REAL(wp)                       ::   zmax3
      INTEGER  ::   jlay, ji, jj, jk, jl   ! dummy loop indices
      !!----------------------------------------------------------------------
      jlay = SIZE( pice , 3 )   ! size of input arrays
      ! basic version: get the max of epsi20 + 9 neighbours
!!$      DO jl = 1, jpl
!!$         DO jk = 1, jlay
!!$            DO_2D( 0, 0, 0, 0 )
!!$               pmax(ji,jj,jk,jl) = MAX( epsi20, pice(ji-1,jj-1,jk,jl), pice(ji,jj-1,jk,jl), pice(ji+1,jj-1,jk,jl),   &
!!$                  &                             pice(ji-1,jj  ,jk,jl), pice(ji,jj  ,jk,jl), pice(ji+1,jj  ,jk,jl),   &
!!$                  &                             pice(ji-1,jj+1,jk,jl), pice(ji,jj+1,jk,jl), pice(ji+1,jj+1,jk,jl) )
!!$            END_2D
!!$         END DO
!!$      END DO
      ! optimized version : does a little bit more than 2 max of epsi20 + 3 neighbours
      DO jl = 1, jpl
         DO jk = 1, jlay
            DO ji = Nis0, Nie0
               zmax1(ji) = MAX( epsi20, pice(ji,Njs0-1,jk,jl), pice(ji-1,Njs0-1,jk,jl), pice(ji+1,Njs0-1,jk,jl) )
               zmax2(ji) = MAX( epsi20, pice(ji,Njs0  ,jk,jl), pice(ji-1,Njs0  ,jk,jl), pice(ji+1,Njs0  ,jk,jl) )
            END DO
            DO_2D( 0, 0, 0, 0 )
               zmax3 = MAX( epsi20, pice(ji,jj+1,jk,jl), pice(ji-1,jj+1,jk,jl), pice(ji+1,jj+1,jk,jl) )
               pmax(ji,jj,jk,jl) = MAX( epsi20, zmax1(ji), zmax2(ji), zmax3 )
               zmax1(ji) = zmax2(ji)
               zmax2(ji) = zmax3
            END_2D
         END DO
      END DO
   END SUBROUTINE icemax4D

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

   !!======================================================================
END MODULE icedyn_adv_pra