source: tags/LIM1D_v3.20/SOURCES/source_3.20/ice_phy_ini.f @ 6

      SUBROUTINE ice_phy_ini(kideb,kiut, nlay_s, nlay_i)

      ! BIO-LIM1D, 2008 (M. Vancoppenolle)
      ! update 2014, v3.04

      ! This routine initializes ice variables
      
!------------------------------------------------------------------------------!

      INCLUDE 'type.com'
      INCLUDE 'para.com'
      INCLUDE 'const.com'
      INCLUDE 'ice.com'
      INCLUDE 'thermo.com'
      INCLUDE 'bio.com'

      CHARACTER(len=10) :: 
     &   filenc='init.nc'

      REAL(4)        zini(1)  ! forcing field dummy array
      REAL(4)        zzi(maxnlay), zti(maxnlay), zsi(maxnlay) ! forcing field dummy array
      REAL(4)        zdsi(maxnlay), zno3(maxnlay), zpo4(maxnlay)
      REAL(4)        zdaf(maxnlay)

      REAL(8), DIMENSION (maxnlay) :: 
     &   zsh_i0              ,   !: old ice salt content (ppt.m-2)
     &   zsh_i1                  !: new ice salt content (ppt.m-2)

      REAL(8), DIMENSION (maxnlay) :: 
     &   zqh_i0              ,   !: old ice heat content (J.m-2)
     &   zqh_i1                  !: new ice heat content (J.m-2)

      REAL(8), DIMENSION (maxnlay+2) ::
     &   zthick0                 !: thickness of old layers

      WRITE(numout,*) ' * ice_phy_ini : '
      WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
      WRITE(numout,*) 

      zeps = 1.0e-10
      ji = 1

!
!------------------------------------------------------------------------------|
!  1) Case of a NETCDF file
!------------------------------------------------------------------------------|
!
      !-----------
      ! Open file
      !-----------
      CALL CF_OPEN  (filenc,id) ! open forcing file

      !-----------
      ! Read data
      !-----------
      zn = 0.0
      CALL CF_READ1D ( filenc, 'nlay', 1, 1, zini )
      zn = REAL(zini(1))

      zj_d = 0.0
      CALL CF_READ1D ( filenc, 'j_d', 1, 1, zini )
      zj_d = REAL(zini(1))

      CALL CF_READ1D ( filenc, 'h_s', 1, 1, zini)
      ht_s_b(ji) = REAL(zini(1))
       

      CALL CF_READ1D ( filenc, 'h_i', 1, 1, zini)
      ht_i_b(ji) = REAL(zini(1))

      CALL CF_READ1D ( filenc, 't_s', 1, 1, zini)
      t_s_b(ji,1) = REAL(zini(1))

      CALL CF_READ1D ( filenc, 't_su', 1, 1, zini)
      t_su_b(ji) = REAL(zini(1))

      DO layer = 1, INT(zn)
         CALL CF_READ1D ( filenc, 'z_i', layer, 1, zini)
         zzi(layer) = zini(1)
         CALL CF_READ1D ( filenc, 't_i', layer, 1, zini)
         zti(layer) = zini(1)
         CALL CF_READ1D ( filenc, 's_i', layer, 1, zini)
         zsi(layer) = zini(1)
      END DO

      CALL CF_CLOSE (filenc) ! close nc file

      t_bo_b(ji)   =  273.15-tmut*oce_sal 

      !----------------------
      ! Interpolate profiles
      !----------------------

            !------
            ! Snow
            !------
            CALL ice_phy_grid(1,1,nlay_s,ht_s_b(ji), .FALSE., "sno" )   ! compute the physical grid

            !-----
            ! Ice
            !-----
            CALL ice_phy_grid(1,1,nlay_i,ht_i_b(ji), .FALSE., "ice" )   ! compute the physical grid

            ! grid indexes for redistribution
            ntop0  = 1 
            nbot0  = INT(zn)
            ntop1  = 1
            nbot1  = nlay_i

            !--- Salinity
            zm0(0) = 0.0         ! layer interfaces cotes
            DO layer = 1, nbot0-1
               zm0(layer) = REAL( zzi(layer) + zzi(layer+1) ) / 2.
            END DO
            zm0(nbot0) = ht_i_b(ji)

            DO layer = 1, nbot0  !layer thickness and salt content
               zthick0(layer) = zm0(layer)-zm0(layer-1)
               zsh_i0(layer) = rhog*zthick0(layer)*REAL(zsi(layer)) !mass of salt
            END DO

            CALL ice_phy_relay( nbot0 , nbot1 , ntop0 , ntop1 , 
     &                       zthick0, deltaz_i_phy, zsh_i0 , zsh_i1 )

            DO layer = 1, nlay_i
               s_i_b(ji,layer) = zsh_i1(layer) / 
     &         (rhog*deltaz_i_phy(layer))
            END DO

            DO layer = 1, nbot0  ! heat content
               ztmelts    =   - tmut * REAL(zsi(layer)) + tpw 
               zq         = rhog * ( cpg * ( ztmelts -REAL(zti(layer)) )
     &                      + lfus*( 1.0 - (ztmelts-tpw) / 
     &                        MIN( (REAL( zti(layer) ) -tpw),-zeps) )  
     &                      - cpw      * ( ztmelts-tpw  ) ) 
               zqh_i0(layer) = zq * zthick0(layer)
           END DO
           CALL ice_phy_relay( nbot0 , nbot1 , ntop0 , ntop1 , 
     &                       zthick0, deltaz_i_phy, zqh_i0 , zqh_i1 )

           DO layer = 1, nlay_i
              zq = zqh_i1(layer) / MAX( deltaz_i_phy(layer) , zeps )
              ztmelts = -tmut*s_i_b(ji,layer) + tpw
              aaa = cpg
              bbb = (cpw-cpg)*(ztmelts-tpw) + zq / rhog - lfus
              ccc = lfus * (ztmelts-tpw)
              discrim = SQRT( bbb*bbb - 4.0*aaa*ccc )
              t_i_b(ji,layer) = tpw + (- bbb - discrim) / ( 2.0*aaa )
           END DO

!
!------------------------------------------------------------------------------|
! x) Write in output file
!------------------------------------------------------------------------------|

            WRITE(numout,*) ' --- Initial values --- '
            WRITE(numout,*) ' '
            WRITE(numout,*) ' Surface temperature t_su : ', t_su_b(ji)
            WRITE(numout,*) ' Basal   temperature t_bo : ', t_bo_b(ji)
            WRITE(numout,*) ' Snow temperature    t_s  : ', t_s_b(ji,1)
            WRITE(numout,*) ' Ice temperatures    t_i  : ', 
     &                      ( t_i_b(ji,jk), jk = 1, nlay_i )
            WRITE(numout,*) ' Ice salinities      s_i  : ', 
     &                      ( s_i_b(ji,jk), jk = 1, nlay_i )
            WRITE(numout,*) ' Ice thickness       ht_i_b : ', ht_i_b(ji)
            WRITE(numout,*) ' Snow depth          ht_s_b : ', ht_s_b(ji)
            WRITE(numout,*) ' '

!------------------------------------------------------------------------------|
!- end of ice_phy_ini
!
      RETURN
      END