New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 15334 for NEMO/trunk/src/ICE/icethd.F90 – NEMO

Ignore:
Timestamp:
2021-10-05T23:18:34+02:00 (14 months ago)
Author:
clem
Message:

Some harmless reorganization of SI3: 1) extract the parts where mpi com were needed from inside thermo. 2) code an optional upstream scheme inside rheology to calculate P as the sub time step level. 3) prepare the albedo to scheme to recieve an aditional namelist parameter (pivotal ice thickness)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/ICE/icethd.F90

    r14997 r15334  
    2020   USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, sprecip, ln_cpl 
    2121   USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, & 
    22       &                 qml_ice, qcn_ice, qtr_ice_top 
     22      &                 qml_ice, qcn_ice, qtr_ice_top, utau_ice, vtau_ice 
    2323   USE ice1D          ! sea-ice: thermodynamics variables 
    2424   USE icethd_zdf     ! sea-ice: vertical heat diffusion 
     
    9797      REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io, zfric, zvel   ! ice-ocean velocity (m/s) and frictional velocity (m2/s2) 
    9898      ! 
     99      ! for collection thickness 
     100      INTEGER  ::   iter             !   -       - 
     101      REAL(wp) ::   zvfrx, zvgx, ztaux, zf                     !   -      - 
     102      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, ztwogp   !   -      - 
     103      REAL(wp), PARAMETER ::   zcai = 1.4e-3_wp               ! ice-air drag (clem: should be dependent on coupling/forcing used) 
     104      REAL(wp), PARAMETER ::   zhicrit = 0.04_wp                                             ! frazil ice thickness 
     105      REAL(wp), PARAMETER ::   zsqcd   = 1.0_wp / SQRT( 1.3_wp * zcai )                      ! 1/SQRT(airdensity*drag) 
     106      REAL(wp), PARAMETER ::   zgamafr = 0.03_wp 
    99107      !!------------------------------------------------------------------- 
    100108      ! controls 
     
    128136               &                         ( v_ice(ji,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji,jj-1) + v_ice(ji,jj) ) ) 
    129137         END_2D 
    130       ELSE      !  if no ice dynamics => transmit directly the atmospheric stress to the ocean 
     138      ELSE      !  if no ice dynamics => transfer directly the atmospheric stress to the ocean 
    131139         DO_2D( 0, 0, 0, 0 ) 
    132140            zfric(ji,jj) = r1_rho0 * SQRT( 0.5_wp *  & 
     
    219227      ENDIF 
    220228 
     229 
     230      !---------------------------------------------------------! 
     231      ! Collection thickness of ice formed in leads and polynyas 
     232      !---------------------------------------------------------!     
     233      ! ht_i_new is the thickness of new ice formed in open water 
     234      ! ht_i_new can be either prescribed (ln_frazil=F) or computed (ln_frazil=T) 
     235      ! Frazil ice forms in open water, is transported by wind, accumulates at the edge of the consolidated ice edge 
     236      ! where it forms aggregates of a specific thickness called collection thickness. 
     237      ! 
     238      fraz_frac(:,:) = 0._wp 
     239      ! 
     240      ! Default new ice thickness 
     241      WHERE( qlead(:,:) < 0._wp ) ! cooling 
     242         ht_i_new(:,:) = rn_hinew 
     243      ELSEWHERE 
     244         ht_i_new(:,:) = 0._wp 
     245      END WHERE 
     246 
     247      IF( ln_frazil ) THEN 
     248         ztwogp  = 2._wp * rho0 / ( grav * 0.3_wp * ( rho0 - rhoi ) )  ! reduced grav 
     249         ! 
     250         DO_2D( 0, 0, 0, 0 ) 
     251            IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling 
     252               ! -- Wind stress -- ! 
     253               ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   & 
     254                  &          +   utau_ice(ji  ,jj  ) * umask(ji  ,jj  ,1) ) * 0.5_wp 
     255               ztauy         = ( vtau_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)   & 
     256                  &          +   vtau_ice(ji  ,jj  ) * vmask(ji  ,jj  ,1) ) * 0.5_wp 
     257               ! Square root of wind stress 
     258               ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 
     259 
     260               ! -- Frazil ice velocity -- ! 
     261               rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 
     262               zvfrx   = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 
     263               zvfry   = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 
     264 
     265               ! -- Pack ice velocity -- ! 
     266               zvgx    = ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 
     267               zvgy    = ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 
     268 
     269               ! -- Relative frazil/pack ice velocity -- ! 
     270               rswitch      = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
     271               zvrel2       = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   & 
     272                  &               + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15_wp * 0.15_wp ) * rswitch 
     273 
     274               !!clem 
     275               fraz_frac(ji,jj) = rswitch * ( TANH( rn_Cfraz * ( SQRT(zvrel2) - rn_vfraz ) ) + 1._wp ) * 0.5_wp * rn_maxfraz 
     276               !!clem 
     277                
     278               ! -- new ice thickness (iterative loop) -- ! 
     279               ht_i_new(ji,jj) = zhicrit +   ( zhicrit + 0.1_wp )    & 
     280                  &                      / ( ( zhicrit + 0.1_wp ) * ( zhicrit + 0.1_wp ) -  zhicrit * zhicrit ) * ztwogp * zvrel2 
     281               iter = 1 
     282               DO WHILE ( iter < 20 )  
     283                  zf  = ( ht_i_new(ji,jj) - zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj) - zhicrit * zhicrit ) -   & 
     284                     &    ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2 
     285                  zfp = ( ht_i_new(ji,jj) - zhicrit ) * ( 3.0_wp * ht_i_new(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2 
     286 
     287                  ht_i_new(ji,jj) = ht_i_new(ji,jj) - zf / MAX( zfp, epsi20 ) 
     288                  iter = iter + 1 
     289               END DO 
     290               ! 
     291               ! bound ht_i_new (though I don't see why it should be necessary) 
     292               ht_i_new(ji,jj) = MAX( 0.01_wp, MIN( ht_i_new(ji,jj), hi_max(jpl) ) ) 
     293               ! 
     294            ELSE 
     295               ht_i_new(ji,jj) = 0._wp 
     296            ENDIF 
     297            ! 
     298         END_2D 
     299         !  
     300         CALL lbc_lnk( 'icethd', fraz_frac, 'T', 1.0_wp, ht_i_new, 'T', 1.0_wp  ) 
     301 
     302      ENDIF 
     303 
    221304      !-------------------------------------------------------------------------------------------! 
    222305      ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories 
     
    268351      ! 
    269352      IF ( ln_pnd .AND. ln_icedH ) & 
    270          &                    CALL ice_thd_pnd                      ! --- Melt ponds 
     353         &                    CALL ice_thd_pnd                      ! --- Melt ponds --- ! 
    271354      ! 
    272355      IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- ! 
     
    276359                              CALL ice_cor( kt , 2 )                ! --- Corrections --- ! 
    277360      ! 
    278       oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice              ! ice natural aging incrementation 
     361      oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice              ! --- Ice natural aging incrementation 
     362      ! 
     363      DO_2D( 0, 0, 0, 0 )                                           ! --- Ice velocity corrections 
     364         IF( at_i(ji,jj) == 0._wp ) THEN   ! if ice has melted 
     365            IF( at_i(ji+1,jj) == 0._wp )   u_ice(ji  ,jj) = 0._wp   ! right side 
     366            IF( at_i(ji-1,jj) == 0._wp )   u_ice(ji-1,jj) = 0._wp   ! left side 
     367            IF( at_i(ji,jj+1) == 0._wp )   v_ice(ji,jj  ) = 0._wp   ! upper side 
     368            IF( at_i(ji,jj-1) == 0._wp )   v_ice(ji,jj-1) = 0._wp   ! bottom side 
     369         ENDIF 
     370      END_2D 
     371      CALL lbc_lnk( 'icecor', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) 
    279372      ! 
    280373      ! convergence tests 
Note: See TracChangeset for help on using the changeset viewer.