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 – NEMO

Changeset 15334


Ignore:
Timestamp:
2021-10-05T23:18:34+02:00 (3 years 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)

Location:
NEMO/trunk/src/ICE
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/src/ICE/ice.F90

    r14103 r15334  
    147147   ! 
    148148   !                                     !!** ice-ridging/rafting namelist (namdyn_rdgrft) ** 
     149   LOGICAL,  PUBLIC ::   ln_str_H79       !: ice strength parameterization (Hibler79) (may be used in rheology) 
    149150   REAL(wp), PUBLIC ::   rn_crhg          !: determines changes in ice strength (also used for landfast param) 
     151   REAL(wp), PUBLIC ::   rn_pstar         !: determines ice strength, Hibler JPO79 (may be used in rheology) 
    150152   ! 
    151153   !                                     !!** ice-rheology namelist (namdyn_rhg) ** 
    … …  
    251253   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_oce,v_oce     !: surface ocean velocity used in ice dynamics 
    252254   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ht_i_new        !: ice collection thickness accreted in leads 
     255   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fraz_frac       !: fraction of frazil ice accreted at the ice bottom 
    253256   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   strength        !: ice strength 
    254257   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   stress1_i, stress2_i, stress12_i   !: 1st, 2nd & diagonal stress tensor element 
    … …  
    453456 
    454457      ii = 1 
    455       ALLOCATE( u_oce    (jpi,jpj) , v_oce    (jpi,jpj) , ht_i_new  (jpi,jpj) , strength(jpi,jpj) ,  & 
    456          &      stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,                      & 
    457          &      delta_i  (jpi,jpj) , divu_i   (jpi,jpj) , shear_i   (jpi,jpj) ,                      & 
    458          &      aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , rdg_conv  (jpi,jpj) , STAT=ierr(ii) ) 
     458      ALLOCATE( u_oce    (jpi,jpj) , v_oce    (jpi,jpj) , ht_i_new (jpi,jpj) , fraz_frac (jpi,jpj) ,  & 
     459         &      strength (jpi,jpj) , stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,  & 
     460         &      delta_i  (jpi,jpj) , divu_i   (jpi,jpj) , shear_i  (jpi,jpj) ,                        & 
     461         &      aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , rdg_conv (jpi,jpj) , STAT=ierr(ii) ) 
    459462 
    460463      ii = ii + 1 

  • NEMO/trunk/src/ICE/icealb.F90

    r14997 r15334  
    3030   PUBLIC   ice_alb        ! called in icesbc.F90 and iceupdate.F90 
    3131 
    32    REAL(wp), PUBLIC, PARAMETER ::   rn_alb_oce = 0.066   !: ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001) 
     32   REAL(wp), PUBLIC, PARAMETER ::   rn_alb_oce = 0.066_wp   !: ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001) 
    3333   ! 
    3434   !                             !!* albedo namelist (namalb) 
    … …  
    111111      REAL(wp) ::   zalb_snw, zafrac_snw      ! snow-covered sea ice albedo & relative snow fraction 
    112112      REAL(wp) ::   zalb_cs, zalb_os          ! albedo of ice under clear/overcast sky 
     113      !! clem 
     114      REAL(wp), PARAMETER ::   zhi_albcst = 1.5_wp ! pivotal thickness (should be in the namelist) 
    113115      !!--------------------------------------------------------------------- 
    114116      ! 
    115117      IF( ln_timing )   CALL timing_start('icealb') 
    116118      ! 
    117       z1_href_pnd = 1. / 0.05 
    118       z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) )  
    119       z1_c2 = 1. / 0.05 
    120       z1_c3 = 1. / 0.02 
    121       z1_c4 = 1. / 0.03 
     119      z1_href_pnd = 1._wp / 0.05_wp 
     120      z1_c1 = 1._wp / ( LOG(zhi_albcst) - LOG(0.05_wp) )  
     121      z1_c2 = 1._wp / 0.05_wp 
     122      z1_c3 = 1._wp / 0.02_wp 
     123      z1_c4 = 1._wp / 0.03_wp 
    122124      ! 
    123125      CALL ice_var_snwfra( ph_snw, za_s_fra )   ! calculate ice fraction covered by snow 
    … …  
    148150            ENDIF 
    149151            !                       !--- Bare ice albedo (for hi < 150cm) 
    150             IF( 0.05 < ph_ice(ji,jj,jl) .AND. ph_ice(ji,jj,jl) <= 1.5 ) THEN      ! 5cm < hi < 150cm 
    151                zalb_ice = zalb_ice    + ( 0.18 - zalb_ice   ) * z1_c1 * ( LOG(1.5) - LOG(ph_ice(ji,jj,jl)) ) 
    152             ELSEIF( ph_ice(ji,jj,jl) <= 0.05 ) THEN                               ! 0cm < hi < 5cm 
    153                zalb_ice = rn_alb_oce  + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice(ji,jj,jl) 
     152            IF( 0.05 < ph_ice(ji,jj,jl) .AND. ph_ice(ji,jj,jl) <= zhi_albcst ) THEN      ! 5cm < hi < 150cm 
     153               zalb_ice = zalb_ice    + ( 0.18_wp - zalb_ice   ) * z1_c1 * ( LOG(zhi_albcst) - LOG(ph_ice(ji,jj,jl)) ) 
     154            ELSEIF( ph_ice(ji,jj,jl) <= 0.05_wp ) THEN                               ! 0cm < hi < 5cm 
     155               zalb_ice = rn_alb_oce  + ( 0.18_wp - rn_alb_oce ) * z1_c2 * ph_ice(ji,jj,jl) 
    154156            ENDIF 
    155157            ! 
    … …  
    166168            zalb_os = ( zafrac_snw * zalb_snw + zafrac_pnd * zalb_pnd + zafrac_ice * zalb_ice ) * tmask(ji,jj,1) 
    167169            ! 
    168             zalb_cs = zalb_os - ( - 0.1010 * zalb_os * zalb_os  & 
    169                &                  + 0.1933 * zalb_os - 0.0148 ) * tmask(ji,jj,1) 
     170            zalb_cs = zalb_os - ( - 0.1010_wp * zalb_os * zalb_os  & 
     171               &                  + 0.1933_wp * zalb_os - 0.0148_wp ) * tmask(ji,jj,1) 
    170172            ! 
    171173            ! albedo depends on cloud fraction because of non-linear spectral effects 

  • NEMO/trunk/src/ICE/icecor.F90

    r14997 r15334  
    5353      INTEGER, INTENT(in) ::   kn    ! 1 = after dyn ; 2 = after thermo 
    5454      ! 
    55       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
     55      INTEGER  ::   ji, jj, jl       ! dummy loop indices 
    5656      REAL(wp) ::   zsal, zzc 
    5757      !!---------------------------------------------------------------------- 
    … …  
    9999         END DO 
    100100      ENDIF 
    101  
     101      ! 
    102102      IF( kn /= 0 ) THEN   ! no zapsmall if kn=0 (for bdy for instance) because we do not want ice-ocean exchanges (wfx,sfx,hfx) 
    103103         !                                                              otherwise conservation diags will fail 
    … …  
    105105         CALL ice_var_zapsmall      !  Zap small values                                  ! 
    106106         !                          !----------------------------------------------------- 
    107       ENDIF 
    108       !                             !----------------------------------------------------- 
    109       IF( kn == 2 ) THEN            !  Ice drift case: Corrections to avoid wrong values ! 
    110          DO_2D( 0, 0, 0, 0 )        !----------------------------------------------------- 
    111             IF ( at_i(ji,jj) == 0._wp ) THEN    ! what to do if there is no ice 
    112                IF ( at_i(ji+1,jj) == 0._wp )   u_ice(ji  ,jj) = 0._wp   ! right side 
    113                IF ( at_i(ji-1,jj) == 0._wp )   u_ice(ji-1,jj) = 0._wp   ! left side 
    114                IF ( at_i(ji,jj+1) == 0._wp )   v_ice(ji,jj  ) = 0._wp   ! upper side 
    115                IF ( at_i(ji,jj-1) == 0._wp )   v_ice(ji,jj-1) = 0._wp   ! bottom side 
    116             ENDIF 
    117          END_2D 
    118          CALL lbc_lnk( 'icecor', u_ice, 'U', -1.0_wp, v_ice, 'V', -1.0_wp ) 
    119107      ENDIF 
    120108      ! 

  • NEMO/trunk/src/ICE/icedyn_rdgrft.F90

    r14997 r15334  
    5757   ! 
    5858   ! ** namelist (namdyn_rdgrft) ** 
    59    LOGICAL  ::   ln_str_H79       ! ice strength parameterization (Hibler79) 
    60    REAL(wp) ::   rn_pstar         ! determines ice strength, Hibler JPO79 
    6159   REAL(wp) ::   rn_csrdg         ! fraction of shearing energy contributing to ridging 
    6260   LOGICAL  ::   ln_partf_lin     ! participation function linear (Thorndike et al. (1975)) 

  • NEMO/trunk/src/ICE/icedyn_rhg_evp.F90

    r15266 r15334  
    5252   INTEGER ::   nvarid   ! netcdf variable id 
    5353   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   fimask   ! mask at F points for the ice 
    54  
     54    
    5555   !! * Substitutions 
    5656#  include "do_loop_substitute.h90" 
    … …  
    180180      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s) 
    181181      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s) 
     182      !! -- advect fields at the rheology time step for the calculation of strength 
     183      LOGICAL  ::   ll_advups = .FALSE. 
     184      REAL(wp) ::   zdt_ups 
     185      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   ztmp 
     186      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   za_i_ups, zv_i_ups   ! tracers advected upstream 
    182187      !!------------------------------------------------------------------- 
    183188 
    … …  
    234239         z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
    235240      ELSE 
    236          zdtevp   = rdt_ice 
     241         zdtevp   = rDt_ice 
    237242         ! zalpha parameters set later on adaptatively 
    238243      ENDIF 
    … …  
    383388            zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 
    384389 
    385          END_2D 
    386          CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1.0_wp ) 
    387  
    388          ! P/delta at T points 
    389          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     390            ! P/delta at T points 
    390391            zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl ) 
    391          END_2D 
    392  
     392 
     393         END_2D 
     394         CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1.0_wp, zp_delt, 'T', 1.0_wp ) 
     395 
     396         ! 
    393397         DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )   ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2 
    394398 
    … …  
    686690         IF( nn_rhg_chkcvg == 2 )   CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice, zmsk15 ) 
    687691         ! 
     692         ! 
     693         ! --- change strength according to advected a_i and v_i (upstream for now) --- ! 
     694         IF( ll_advups .AND. ln_str_H79 ) THEN 
     695            ! 
     696            IF( jter == 1 ) THEN                               ! init 
     697               ALLOCATE( za_i_ups(jpi,jpj,jpl), zv_i_ups(jpi,jpj,jpl) ) 
     698               ALLOCATE( ztmp(jpi,jpj) ) 
     699               zdt_ups = rDt_ice / REAL( nn_nevp ) 
     700               za_i_ups(:,:,:) = a_i(:,:,:) 
     701               zv_i_ups(:,:,:) = v_i(:,:,:) 
     702            ELSE 
     703               CALL lbc_lnk( 'icedyn_rhg_evp', za_i_ups, 'T', 1.0_wp, zv_i_ups, 'T', 1.0_wp )                
     704            ENDIF 
     705            ! 
     706            CALL rhg_upstream( zdt_ups, u_ice, v_ice, za_i_ups )   ! upstream advection: a_i 
     707            CALL rhg_upstream( zdt_ups, u_ice, v_ice, zv_i_ups )   ! upstream advection: v_i 
     708            ! 
     709            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )    ! strength 
     710               strength(ji,jj) = rn_pstar * SUM( zv_i_ups(ji,jj,:) ) * EXP( -rn_crhg * ( 1._wp - SUM( za_i_ups(ji,jj,:) ) ) ) 
     711            END_2D 
     712            IF( nn_hls == 1 )  CALL lbc_lnk( 'icedyn_rhg_evp', strength, 'T', 1.0_wp ) 
     713            ! 
     714            DO_2D( 0, 0, 0, 0 )                                ! strength smoothing 
     715               IF( SUM( za_i_ups(ji,jj,:) ) > 0._wp ) THEN 
     716                  ztmp(ji,jj) = ( 4._wp * strength(ji,jj) + strength(ji-1,jj  ) + strength(ji+1,jj  ) & 
     717                     &                                    + strength(ji  ,jj-1) + strength(ji  ,jj+1) & 
     718                     &          ) / ( 4._wp + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) ) 
     719               ELSE 
     720                  ztmp(ji,jj) = 0._wp 
     721               ENDIF 
     722            END_2D 
     723            DO_2D( 0, 0, 0, 0 ) 
     724               strength(ji,jj) = ztmp(ji,jj) 
     725            END_2D 
     726            ! 
     727            IF( jter == nn_nevp ) THEN 
     728               DEALLOCATE( za_i_ups, zv_i_ups ) 
     729               DEALLOCATE( ztmp ) 
     730            ENDIF 
     731         ENDIF 
    688732         !                                                ! ==================== ! 
    689733      END DO                                              !  end loop over jter  ! 
    690734      !                                                   ! ==================== ! 
    691735      IF( ln_aEVP )   CALL iom_put( 'beta_evp' , zbeta ) 
     736      ! 
     737      IF( ll_advups .AND. ln_str_H79 )   CALL lbc_lnk( 'icedyn_rhg_evp', strength, 'T', 1.0_wp ) 
    692738      ! 
    693739      !------------------------------------------------------------------------------! 
    … …  
    10231069   END SUBROUTINE rhg_evp_rst 
    10241070 
     1071   SUBROUTINE rhg_upstream( pdt, pu, pv, pt ) 
     1072      !!--------------------------------------------------------------------- 
     1073      !!                    ***  ROUTINE rhg_upstream  *** 
     1074      !! 
     1075      !! **  Purpose :   compute the upstream fluxes and upstream guess of tracer 
     1076      !!---------------------------------------------------------------------- 
     1077      REAL(wp)                   , INTENT(in   ) ::   pdt              ! tracer time-step 
     1078      REAL(wp), DIMENSION(:,:  ) , INTENT(in   ) ::   pu, pv           ! 2 ice velocity components 
     1079      REAL(wp), DIMENSION(:,:,:) , INTENT(inout) ::   pt               ! tracer fields 
     1080      ! 
     1081      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
     1082      REAL(wp) ::   ztra          ! local scalar 
     1083      REAL(wp), DIMENSION(jpi,jpj) ::   zfu_ups, zfv_ups   ! upstream fluxes 
     1084      !!---------------------------------------------------------------------- 
     1085      DO jl = 1, jpl 
     1086         DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 
     1087            zfu_ups(ji,jj) = MAX( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji  ,jj  ,jl) + & 
     1088               &             MIN( pu(ji,jj)*e2u(ji,jj), 0._wp ) * pt(ji+1,jj  ,jl) 
     1089            zfv_ups(ji,jj) = MAX( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji  ,jj  ,jl) + & 
     1090               &             MIN( pv(ji,jj)*e1v(ji,jj), 0._wp ) * pt(ji  ,jj+1,jl) 
     1091         END_2D 
     1092         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 
     1093            ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) + zfv_ups(ji,jj) - zfv_ups(ji,jj-1) ) 
     1094            ! 
     1095            pt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 
     1096         END_2D 
     1097      END DO 
     1098      ! 
     1099   END SUBROUTINE rhg_upstream 
    10251100 
    10261101#else 

  • 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 

  • NEMO/trunk/src/ICE/icethd_do.F90

    r14997 r15334  
    1717   USE phycst         ! physical constants 
    1818   USE sbc_oce , ONLY : sss_m 
    19    USE sbc_ice , ONLY : utau_ice, vtau_ice 
    2019   USE ice1D          ! sea-ice: thermodynamics variables 
    2120   USE ice            ! sea-ice: variables 
    … …  
    3837 
    3938   !                          !!** namelist (namthd_do) ** 
    40    REAL(wp) ::   rn_hinew      ! thickness for new ice formation (m) 
    41    LOGICAL ::   ln_frazil     ! use of frazil ice collection as function of wind (T) or not (F) 
    42    REAL(wp) ::   rn_maxfraz    ! maximum portion of frazil ice collecting at the ice bottom 
    43    REAL(wp) ::   rn_vfraz      ! threshold drift speed for collection of bottom frazil ice 
    44    REAL(wp) ::   rn_Cfraz      ! squeezing coefficient for collection of bottom frazil ice 
     39   REAL(wp), PUBLIC ::   rn_hinew      ! thickness for new ice formation (m) 
     40   LOGICAL , PUBLIC ::   ln_frazil     ! use of frazil ice collection as function of wind (T) or not (F) 
     41   REAL(wp), PUBLIC ::   rn_maxfraz    ! maximum portion of frazil ice collecting at the ice bottom 
     42   REAL(wp), PUBLIC ::   rn_vfraz      ! threshold drift speed for collection of bottom frazil ice 
     43   REAL(wp), PUBLIC ::   rn_Cfraz      ! squeezing coefficient for collection of bottom frazil ice 
    4544 
    4645   !! * Substitutions 
    … …  
    7877      !!------------------------------------------------------------------------ 
    7978      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    80       INTEGER  ::   iter             !   -       - 
    81       REAL(wp) ::   ztmelts, zfrazb, zweight, zde                               ! local scalars 
    82       REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf                     !   -      - 
    83       REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      - 
    84       ! 
     79      ! 
     80      REAL(wp) ::   ztmelts 
     81      REAL(wp) ::   zdE 
    8582      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean) 
    8683      REAL(wp) ::   zEi          ! sea ice specific enthalpy (J/kg) 
    … …  
    10299      REAL(wp), DIMENSION(jpij) ::   zda_res     ! residual area in case of excessive heat budget 
    103100      REAL(wp), DIMENSION(jpij) ::   zv_frazb    ! accretion of frazil ice at the ice bottom 
    104       REAL(wp), DIMENSION(jpij) ::   zvrel_1d    ! relative ice / frazil velocity (1D vector) 
     101      REAL(wp), DIMENSION(jpij) ::   zfraz_frac_1d ! relative ice / frazil velocity (1D vector) 
    105102      ! 
    106103      REAL(wp), DIMENSION(jpij,jpl) ::   zv_b    ! old volume of ice in category jl 
    … …  
    109106      REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_2d !: 1-D version of e_i 
    110107      ! 
    111       REAL(wp), DIMENSION(jpi,jpj) ::   zvrel    ! relative ice / frazil velocity 
    112       ! 
    113       REAL(wp) :: zcai = 1.4e-3_wp               ! ice-air drag (clem: should be dependent on coupling/forcing used) 
    114108      !!-----------------------------------------------------------------------! 
    115109 
    … …  
    119113      at_i(:,:) = SUM( a_i, dim=3 ) 
    120114      !------------------------------------------------------------------------------! 
    121       ! 1) Collection thickness of ice formed in leads and polynyas 
    122       !------------------------------------------------------------------------------!     
    123       ! ht_i_new is the thickness of new ice formed in open water 
    124       ! ht_i_new can be either prescribed (ln_frazil=F) or computed (ln_frazil=T) 
    125       ! Frazil ice forms in open water, is transported by wind 
    126       ! accumulates at the edge of the consolidated ice edge 
    127       ! where it forms aggregates of a specific thickness called 
    128       ! collection thickness. 
    129  
    130       zvrel(:,:) = 0._wp 
    131  
    132       ! Default new ice thickness 
    133       WHERE( qlead(:,:) < 0._wp ) ! cooling 
    134          ht_i_new(:,:) = rn_hinew 
    135       ELSEWHERE 
    136          ht_i_new(:,:) = 0._wp 
    137       END WHERE 
    138  
    139       IF( ln_frazil ) THEN 
    140          ! 
    141          ht_i_new(:,:) = 0._wp 
    142          ! 
    143          ! Physical constants 
    144          zhicrit = 0.04                                          ! frazil ice thickness 
    145          ztwogp  = 2. * rho0 / ( grav * 0.3 * ( rho0 - rhoi ) )  ! reduced grav 
    146          zsqcd   = 1.0 / SQRT( 1.3 * zcai )                      ! 1/SQRT(airdensity*drag) 
    147          zgamafr = 0.03 
    148          ! 
    149          DO_2D( 0, 0, 0, 0 ) 
    150             IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling 
    151                ! -- Wind stress -- ! 
    152                ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   & 
    153                   &          +   utau_ice(ji  ,jj  ) * umask(ji  ,jj  ,1) ) * 0.5_wp 
    154                ztauy         = ( vtau_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)   & 
    155                   &          +   vtau_ice(ji  ,jj  ) * vmask(ji  ,jj  ,1) ) * 0.5_wp 
    156                ! Square root of wind stress 
    157                ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 
    158  
    159                ! -- Frazil ice velocity -- ! 
    160                rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 
    161                zvfrx   = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 
    162                zvfry   = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 
    163  
    164                ! -- Pack ice velocity -- ! 
    165                zvgx    = ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 
    166                zvgy    = ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 
    167  
    168                ! -- Relative frazil/pack ice velocity -- ! 
    169                rswitch      = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
    170                zvrel2       = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   & 
    171                   &               + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) * rswitch 
    172                zvrel(ji,jj) = SQRT( zvrel2 ) 
    173  
    174                ! -- new ice thickness (iterative loop) -- ! 
    175                ht_i_new(ji,jj) = zhicrit +   ( zhicrit + 0.1 )    & 
    176                   &                   / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) -  zhicrit * zhicrit ) * ztwogp * zvrel2 
    177  
    178                iter = 1 
    179                DO WHILE ( iter < 20 )  
    180                   zf  = ( ht_i_new(ji,jj) - zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj) - zhicrit * zhicrit ) -   & 
    181                      &    ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2 
    182                   zfp = ( ht_i_new(ji,jj) - zhicrit ) * ( 3.0 * ht_i_new(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2 
    183  
    184                   ht_i_new(ji,jj) = ht_i_new(ji,jj) - zf / MAX( zfp, epsi20 ) 
    185                   iter = iter + 1 
    186                END DO 
    187                ! 
    188                ! bound ht_i_new (though I don't see why it should be necessary) 
    189                ht_i_new(ji,jj) = MAX( 0.01_wp, MIN( ht_i_new(ji,jj), hi_max(jpl) ) ) 
    190                ! 
    191             ENDIF 
    192             ! 
    193          END_2D 
    194          !  
    195          CALL lbc_lnk( 'icethd_do', zvrel, 'T', 1.0_wp, ht_i_new, 'T', 1.0_wp  ) 
    196  
    197       ENDIF 
    198  
    199       !------------------------------------------------------------------------------! 
    200       ! 2) Compute thickness, salinity, enthalpy, age, area and volume of new ice 
     115      ! 1) Compute thickness, salinity, enthalpy, age, area and volume of new ice 
    201116      !------------------------------------------------------------------------------! 
    202117      ! it occurs if cooling 
    … …  
    223138            END DO 
    224139         END DO 
    225          CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d  (1:npti) , qlead      ) 
    226          CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d   (1:npti) , t_bo       ) 
    227          CALL tab_2d_1d( npti, nptidx(1:npti), sfx_opw_1d(1:npti) , sfx_opw    ) 
    228          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_opw_1d(1:npti) , wfx_opw    ) 
    229          CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , ht_i_new   ) 
    230          CALL tab_2d_1d( npti, nptidx(1:npti), zvrel_1d  (1:npti) , zvrel      ) 
     140         CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d     (1:npti) , qlead      ) 
     141         CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d      (1:npti) , t_bo       ) 
     142         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_opw_1d   (1:npti) , sfx_opw    ) 
     143         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_opw_1d   (1:npti) , wfx_opw    ) 
     144         CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice    (1:npti) , ht_i_new   ) 
     145         CALL tab_2d_1d( npti, nptidx(1:npti), zfraz_frac_1d(1:npti) , fraz_frac  ) 
    231146 
    232147         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d(1:npti) , hfx_thd    ) 
    … …  
    300215         END DO 
    301216          
    302          zv_frazb(1:npti) = 0._wp 
    303          IF( ln_frazil ) THEN 
    304             ! A fraction zfrazb of frazil ice is accreted at the ice bottom 
    305             DO ji = 1, npti 
    306                rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) ) 
    307                zfrazb        = rswitch * ( TANH( rn_Cfraz * ( zvrel_1d(ji) - rn_vfraz ) ) + 1.0 ) * 0.5 * rn_maxfraz 
    308                zv_frazb(ji)  =         zfrazb   * zv_newice(ji) 
    309                zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 
    310             END DO 
    311          END IF 
     217         ! A fraction fraz_frac of frazil ice is accreted at the ice bottom 
     218         DO ji = 1, npti 
     219            rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) ) 
     220            zv_frazb(ji)  =           zfraz_frac_1d(ji) * rswitch   * zv_newice(ji) 
     221            zv_newice(ji) = ( 1._wp - zfraz_frac_1d(ji) * rswitch ) * zv_newice(ji) 
     222         END DO 
    312223          
    313224         ! --- Area of new ice --- ! 
    … …  
    317228 
    318229         !------------------------------------------------------------------------------! 
    319          ! 3) Redistribute new ice area and volume into ice categories                  ! 
     230         ! 2) Redistribute new ice area and volume into ice categories                  ! 
    320231         !------------------------------------------------------------------------------! 
    321232 
Note: See TracChangeset for help on using the changeset viewer.