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 10419 for NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn.F90 – NEMO

Ignore:
Timestamp:
2018-12-19T20:46:30+01:00 (5 years ago)
Author:
smasson
Message:

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10418, see #2133

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn.F90

    r10170 r10419  
    2121   USE icecor         ! sea-ice: corrections 
    2222   USE icevar         ! sea-ice: operations 
     23   USE icectl         ! sea-ice: control prints 
    2324   ! 
    2425   USE in_out_manager ! I/O manager 
     
    3738   INTEGER ::              nice_dyn   ! choice of the type of dynamics 
    3839   !                                        ! associated indices: 
    39    INTEGER, PARAMETER ::   np_dynFULL    = 1   ! full ice dynamics               (rheology + advection + ridging/rafting + correction) 
     40   INTEGER, PARAMETER ::   np_dynALL     = 1   ! full ice dynamics               (rheology + advection + ridging/rafting + correction) 
    4041   INTEGER, PARAMETER ::   np_dynRHGADV  = 2   ! pure dynamics                   (rheology + advection)  
    41    INTEGER, PARAMETER ::   np_dynADV     = 3   ! only advection w prescribed vel.(rn_uvice + advection) 
     42   INTEGER, PARAMETER ::   np_dynADV1D   = 3   ! only advection 1D - test case from Schar & Smolarkiewicz 1996 
     43   INTEGER, PARAMETER ::   np_dynADV2D   = 4   ! only advection 2D w prescribed vel.(rn_uvice + advection) 
    4244   ! 
    4345   ! ** namelist (namdyn) ** 
    44    LOGICAL  ::   ln_dynFULL       ! full ice dynamics               (rheology + advection + ridging/rafting + correction) 
    45    LOGICAL  ::   ln_dynRHGADV     ! no ridge/raft & no corrections  (rheology + advection) 
    46    LOGICAL  ::   ln_dynADV        ! only advection w prescribed vel.(rn_uvice + advection) 
    47    REAL(wp) ::   rn_uice          !    prescribed u-vel (case np_dynADV) 
    48    REAL(wp) ::   rn_vice          !    prescribed v-vel (case np_dynADV) 
     46   LOGICAL  ::   ln_dynALL        ! full ice dynamics                      (rheology + advection + ridging/rafting + correction) 
     47   LOGICAL  ::   ln_dynRHGADV     ! no ridge/raft & no corrections         (rheology + advection) 
     48   LOGICAL  ::   ln_dynADV1D      ! only advection in 1D w ice convergence (test case from Schar & Smolarkiewicz 1996) 
     49   LOGICAL  ::   ln_dynADV2D      ! only advection in 2D w prescribed vel. (rn_uvice + advection) 
     50   REAL(wp) ::   rn_uice          !    prescribed u-vel (case np_dynADV1D & np_dynADV2D) 
     51   REAL(wp) ::   rn_vice          !    prescribed v-vel (case np_dynADV1D & np_dynADV2D) 
    4952    
    5053   !! * Substitutions 
     
    5356   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
    5457   !! $Id$ 
    55    !! Software governed by the CeCILL license (see ./LICENSE) 
     58   !! Software governed by the CeCILL licence     (./LICENSE) 
    5659   !!---------------------------------------------------------------------- 
    5760CONTAINS 
     
    7174      INTEGER, INTENT(in) ::   kt     ! ice time step 
    7275      !! 
    73       INTEGER ::   ji, jj, jl         ! dummy loop indices 
    74       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhmax 
     76      INTEGER  ::   ji, jj, jl        ! dummy loop indices 
     77      REAL(wp) ::   zcoefu, zcoefv 
     78      REAL(wp),              DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max 
     79      REAL(wp), ALLOCATABLE, DIMENSION(:,:)         ::   zdivu_i 
    7580      !!-------------------------------------------------------------------- 
    7681      ! 
    77       IF( ln_timing )   CALL timing_start('icedyn') 
     82      ! controls 
     83      IF( ln_timing    )   CALL timing_start('icedyn')                                                             ! timing 
     84      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    7885      ! 
    7986      IF( kt == nit000 .AND. lwp ) THEN 
     
    8289         WRITE(numout,*)'~~~~~~~' 
    8390      ENDIF 
    84  
    8591      !                       
    86       IF( ln_landfast ) THEN            !-- Landfast ice parameterization: define max bottom friction 
     92      IF( ln_landfast_home ) THEN      !-- Landfast ice parameterization 
    8793         tau_icebfr(:,:) = 0._wp 
    8894         DO jl = 1, jpl 
    89             WHERE( h_i(:,:,jl) > ht_n(:,:) * rn_gamma )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
    90          END DO 
    91          IF( iom_use('tau_icebfr') )   CALL iom_put( 'tau_icebfr', tau_icebfr )   
    92       ENDIF 
    93  
    94       zhmax(:,:,:) = h_i_b(:,:,:)      !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig) 
     95            WHERE( h_i_b(:,:,jl) > ht_n(:,:) * rn_depfra )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 
     96         END DO 
     97      ENDIF 
     98      ! 
     99      !                                !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig) 
    95100      DO jl = 1, jpl 
    96101         DO jj = 2, jpjm1 
     102            DO ji = fs_2, fs_jpim1 
     103               zhip_max(ji,jj,jl) = MAX( epsi20, h_ip_b(ji,jj,jl), h_ip_b(ji+1,jj  ,jl), h_ip_b(ji  ,jj+1,jl), & 
     104                  &                                                h_ip_b(ji-1,jj  ,jl), h_ip_b(ji  ,jj-1,jl), & 
     105                  &                                                h_ip_b(ji+1,jj+1,jl), h_ip_b(ji-1,jj-1,jl), & 
     106                  &                                                h_ip_b(ji+1,jj-1,jl), h_ip_b(ji-1,jj+1,jl) ) 
     107               zhi_max (ji,jj,jl) = MAX( epsi20, h_i_b (ji,jj,jl), h_i_b (ji+1,jj  ,jl), h_i_b (ji  ,jj+1,jl), & 
     108                  &                                                h_i_b (ji-1,jj  ,jl), h_i_b (ji  ,jj-1,jl), & 
     109                  &                                                h_i_b (ji+1,jj+1,jl), h_i_b (ji-1,jj-1,jl), & 
     110                  &                                                h_i_b (ji+1,jj-1,jl), h_i_b (ji-1,jj+1,jl) ) 
     111               zhs_max (ji,jj,jl) = MAX( epsi20, h_s_b (ji,jj,jl), h_s_b (ji+1,jj  ,jl), h_s_b (ji  ,jj+1,jl), & 
     112                  &                                                h_s_b (ji-1,jj  ,jl), h_s_b (ji  ,jj-1,jl), & 
     113                  &                                                h_s_b (ji+1,jj+1,jl), h_s_b (ji-1,jj-1,jl), & 
     114                  &                                                h_s_b (ji+1,jj-1,jl), h_s_b (ji-1,jj+1,jl) ) 
     115            END DO 
     116         END DO 
     117      END DO 
     118      CALL lbc_lnk_multi( 'icedyn', zhi_max(:,:,:), 'T', 1., zhs_max(:,:,:), 'T', 1., zhip_max(:,:,:), 'T', 1. ) 
     119      ! 
     120      ! 
     121      SELECT CASE( nice_dyn )           !-- Set which dynamics is running 
     122 
     123      CASE ( np_dynALL )           !==  all dynamical processes  ==! 
     124         CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology   
     125         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness 
     126         CALL ice_dyn_rdgrft( kt )                                                 ! -- ridging/rafting  
     127         CALL ice_cor       ( kt , 1 )                                             ! -- Corrections 
     128 
     129      CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==! 
     130         CALL ice_dyn_rhg   ( kt )                                                 ! -- rheology   
     131         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness 
     132         CALL Hpiling                                                              ! -- simple pile-up (replaces ridging/rafting) 
     133 
     134      CASE ( np_dynADV1D )         !==  pure advection ==!   (1D) 
     135         ALLOCATE( zdivu_i(jpi,jpj) ) 
     136         ! --- monotonicity test from Schar & Smolarkiewicz 1996 --- ! 
     137         ! CFL = 0.5 at a distance from the bound of 1/6 of the basin length 
     138         ! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s  
     139         DO jj = 1, jpj 
     140            DO ji = 1, jpi 
     141               zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. ) 
     142               zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. ) 
     143               u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1., zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1) 
     144               v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1., zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1) 
     145            END DO 
     146         END DO 
     147         ! --- 
     148         CALL ice_dyn_adv   ( kt )                                       ! -- advection of ice + correction on ice thickness 
     149 
     150         ! diagnostics: divergence at T points  
     151         DO jj = 2, jpjm1 
    97152            DO ji = 2, jpim1 
    98 !!gm use of MAXVAL here is very probably less efficient than expending the 9 values 
    99                zhmax(ji,jj,jl) = MAX( epsi20, MAXVAL( h_i_b(ji-1:ji+1,jj-1:jj+1,jl) ) ) 
    100             END DO 
    101          END DO 
    102       END DO 
    103       CALL lbc_lnk( 'icedyn', zhmax(:,:,:), 'T', 1. ) 
    104       ! 
    105       ! 
    106       SELECT CASE( nice_dyn )           !-- Set which dynamics is running 
    107  
    108       CASE ( np_dynFULL )          !==  all dynamical processes  ==! 
    109          CALL ice_dyn_rhg   ( kt )                            ! -- rheology   
    110          CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhmax )   ! -- advection of ice + correction on ice thickness 
    111          CALL ice_dyn_rdgrft( kt )                            ! -- ridging/rafting  
    112          CALL ice_cor       ( kt , 1 )                        ! -- Corrections 
    113  
    114       CASE ( np_dynRHGADV  )       !==  no ridge/raft & no corrections ==! 
    115          CALL ice_dyn_rhg   ( kt )                            ! -- rheology   
    116          CALL ice_dyn_adv   ( kt )                            ! -- advection of ice 
    117          CALL Hpiling                                         ! -- simple pile-up (replaces ridging/rafting) 
    118  
    119       CASE ( np_dynADV )           !==  pure advection ==!   (prescribed velocities) 
     153               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     154                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 
     155            END DO 
     156         END DO 
     157         CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 
     158         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) ) 
     159 
     160         DEALLOCATE( zdivu_i ) 
     161 
     162      CASE ( np_dynADV2D )         !==  pure advection ==!   (2D w prescribed velocities) 
     163         ALLOCATE( zdivu_i(jpi,jpj) ) 
    120164         u_ice(:,:) = rn_uice * umask(:,:,1) 
    121165         v_ice(:,:) = rn_vice * vmask(:,:,1) 
    122          !!CALL RANDOM_NUMBER(u_ice(:,:)) 
    123          !!CALL RANDOM_NUMBER(v_ice(:,:)) 
    124          CALL ice_dyn_adv   ( kt )                            ! -- advection of ice 
     166         !CALL RANDOM_NUMBER(u_ice(:,:)) ; u_ice(:,:) = u_ice(:,:) * 0.1 + rn_uice * 0.9 * umask(:,:,1) 
     167         !CALL RANDOM_NUMBER(v_ice(:,:)) ; v_ice(:,:) = v_ice(:,:) * 0.1 + rn_vice * 0.9 * vmask(:,:,1) 
     168         ! --- 
     169         CALL ice_dyn_adv   ( kt )   ;   CALL Hbig( zhi_max, zhs_max, zhip_max )   ! -- advection of ice + correction on ice thickness 
     170 
     171         ! diagnostics: divergence at T points  
     172         DO jj = 2, jpjm1 
     173            DO ji = 2, jpim1 
     174               zdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     175                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) ) * r1_e1e2t(ji,jj) 
     176            END DO 
     177         END DO 
     178         CALL lbc_lnk( 'icedyn', zdivu_i, 'T', 1. ) 
     179         IF( iom_use('icediv') )   CALL iom_put( "icediv" , zdivu_i(:,:) ) 
     180 
     181         DEALLOCATE( zdivu_i ) 
    125182 
    126183      END SELECT 
    127       ! 
    128       IF( ln_timing )   CALL timing_stop('icedyn') 
     184       ! 
     185      ! controls 
     186      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
     187      IF( ln_timing    )   CALL timing_stop ('icedyn')                                                             ! timing 
    129188      ! 
    130189   END SUBROUTINE ice_dyn 
    131190 
    132191 
    133    SUBROUTINE Hbig( phmax ) 
     192   SUBROUTINE Hbig( phi_max, phs_max, phip_max ) 
    134193      !!------------------------------------------------------------------- 
    135194      !!                  ***  ROUTINE Hbig  *** 
    136195      !! 
    137196      !! ** Purpose : Thickness correction in case advection scheme creates 
    138       !!              abnormally tick ice 
    139       !! 
    140       !! ** Method  : 1- check whether ice thickness resulting from advection is 
    141       !!                 larger than the surrounding 9-points before advection 
    142       !!                 and reduce it if a) divergence or b) convergence & at_i>0.8 
    143       !!              2- bound ice thickness with hi_max (99m) 
     197      !!              abnormally tick ice or snow 
     198      !! 
     199      !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points 
     200      !!                 (before advection) and reduce it by adapting ice concentration 
     201      !!              2- check whether snow thickness is larger than the surrounding 9-points 
     202      !!                 (before advection) and reduce it by sending the excess in the ocean 
     203      !!              3- check whether snow load deplets the snow-ice interface below sea level$ 
     204      !!                 and reduce it by sending the excess in the ocean 
     205      !!              4- correct pond fraction to avoid a_ip > a_i 
    144206      !! 
    145207      !! ** input   : Max thickness of the surrounding 9-points 
    146208      !!------------------------------------------------------------------- 
    147       REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phmax   ! max ice thick from surrounding 9-pts 
     209      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
    148210      ! 
    149211      INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    150       REAL(wp) ::   zh, zdv 
     212      REAL(wp) ::   zhip, zhi, zhs, zvs_excess, zfra 
    151213      !!------------------------------------------------------------------- 
    152214      ! 
     
    156218         DO jj = 1, jpj 
    157219            DO ji = 1, jpi 
    158                IF ( v_i(ji,jj,jl) > 0._wp ) THEN  !-- bound to hmax 
     220               IF ( v_i(ji,jj,jl) > 0._wp ) THEN 
    159221                  ! 
    160                   zh  = v_i (ji,jj,jl) / a_i(ji,jj,jl) 
    161                   zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl)   
    162                   ! 
    163                   IF ( ( zdv >  0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. & 
    164                      & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN 
    165                      a_i (ji,jj,jl) = v_i(ji,jj,jl) / MIN( phmax(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m) 
     222                  !                               ! -- check h_ip -- ! 
     223                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
     224                  IF( ln_pnd_H12 .AND. v_ip(ji,jj,jl) > 0._wp ) THEN 
     225                     zhip = v_ip(ji,jj,jl) / MAX( epsi20, a_ip(ji,jj,jl) ) 
     226                     IF( zhip > phip_max(ji,jj,jl) .AND. a_ip(ji,jj,jl) < 0.15 ) THEN 
     227                        a_ip(ji,jj,jl) = v_ip(ji,jj,jl) / phip_max(ji,jj,jl) 
     228                     ENDIF 
    166229                  ENDIF 
    167230                  ! 
     231                  !                               ! -- check h_i -- ! 
     232                  ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 
     233                  zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
     234                  IF( zhi > phi_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN 
     235                     a_i(ji,jj,jl) = v_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) )   !-- bound h_i to hi_max (99 m) 
     236                  ENDIF 
     237                  ! 
     238                  !                               ! -- check h_s -- ! 
     239                  ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 
     240                  zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl) 
     241                  IF( v_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. a_i(ji,jj,jl) < 0.15 ) THEN 
     242                     zfra = a_i(ji,jj,jl) * phs_max(ji,jj,jl) / MAX( v_s(ji,jj,jl), epsi20 ) 
     243                     ! 
     244                     wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_s(ji,jj,jl) - a_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * r1_rdtice 
     245                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 
     246                     ! 
     247                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra 
     248                     v_s(ji,jj,jl)          = a_i(ji,jj,jl) * phs_max(ji,jj,jl) 
     249                  ENDIF            
     250                  ! 
     251                  !                               ! -- check snow load -- ! 
     252                  ! if snow load makes snow-ice interface to deplet below the ocean surface => put the snow excess in the ocean 
     253                  !    this correction is crucial because of the call to routine icecor afterwards which imposes a mini of ice thick. (rn_himin) 
     254                  !    this imposed mini can artificially make the snow very thick (if concentration decreases drastically) 
     255                  zvs_excess = MAX( 0._wp, v_s(ji,jj,jl) - v_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 
     256                  IF( zvs_excess > 0._wp ) THEN 
     257                     zfra = zvs_excess / MAX( v_s(ji,jj,jl), epsi20 ) 
     258                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * r1_rdtice 
     259                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( e_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * r1_rdtice ! W.m-2 <0 
     260                     ! 
     261                     e_s(ji,jj,1:nlay_s,jl) = e_s(ji,jj,1:nlay_s,jl) * zfra 
     262                     v_s(ji,jj,jl)          = v_s(ji,jj,jl) - zvs_excess 
     263                  ENDIF 
     264                   
    168265               ENDIF 
    169266            END DO 
     
    215312      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read 
    216313      !! 
    217       NAMELIST/namdyn/ ln_dynFULL, ln_dynRHGADV, ln_dynADV, rn_uice, rn_vice,  & 
    218          &             rn_ishlat  ,                                            & 
    219          &             ln_landfast, rn_gamma , rn_icebfr, rn_lfrelax 
     314      NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice,  & 
     315         &             rn_ishlat ,                                                           & 
     316         &             ln_landfast_L16, ln_landfast_home, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile 
    220317      !!------------------------------------------------------------------- 
    221318      ! 
     
    233330         WRITE(numout,*) '~~~~~~~~~~~~' 
    234331         WRITE(numout,*) '   Namelist namdyn:' 
    235          WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynFULL   = ', ln_dynFULL 
    236          WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV = ', ln_dynRHGADV 
    237          WRITE(numout,*) '      Advection only         (rn_uvice + adv)                 ln_dynADV    = ', ln_dynADV 
    238          WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',', rn_vice,')' 
    239          WRITE(numout,*) '      lateral boundary condition for sea ice dynamics         rn_ishlat    = ', rn_ishlat 
    240          WRITE(numout,*) '      Landfast: param (T or F)                                ln_landfast  = ', ln_landfast 
    241          WRITE(numout,*) '         fraction of ocean depth that ice must reach          rn_gamma     = ', rn_gamma 
    242          WRITE(numout,*) '         maximum bottom stress per unit area of contact       rn_icebfr    = ', rn_icebfr 
    243          WRITE(numout,*) '         relax time scale (s-1) to reach static friction      rn_lfrelax   = ', rn_lfrelax 
     332         WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynALL       = ', ln_dynALL 
     333         WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV    = ', ln_dynRHGADV 
     334         WRITE(numout,*) '      Advection 1D only      (Schar & Smolarkiewicz 1996)     ln_dynADV1D     = ', ln_dynADV1D 
     335         WRITE(numout,*) '      Advection 2D only      (rn_uvice + adv)                 ln_dynADV2D     = ', ln_dynADV2D 
     336         WRITE(numout,*) '         with prescribed velocity given by   (u,v)_ice = (rn_uice,rn_vice)    = (', rn_uice,',', rn_vice,')' 
     337         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics         rn_ishlat       = ', rn_ishlat 
     338         WRITE(numout,*) '      Landfast: param from Lemieux 2016                       ln_landfast_L16 = ', ln_landfast_L16 
     339         WRITE(numout,*) '      Landfast: param from home made                          ln_landfast_home= ', ln_landfast_home 
     340         WRITE(numout,*) '         fraction of ocean depth that ice must reach          rn_depfra       = ', rn_depfra 
     341         WRITE(numout,*) '         maximum bottom stress per unit area of contact       rn_icebfr       = ', rn_icebfr 
     342         WRITE(numout,*) '         relax time scale (s-1) to reach static friction      rn_lfrelax      = ', rn_lfrelax 
     343         WRITE(numout,*) '         isotropic tensile strength                           rn_tensile      = ', rn_tensile 
    244344         WRITE(numout,*) 
    245345      ENDIF 
     
    247347      ioptio = 0  
    248348      !      !--- full dynamics                               (rheology + advection + ridging/rafting + correction) 
    249       IF( ln_dynFULL   ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynFULL      ;   ENDIF 
     349      IF( ln_dynALL    ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynALL       ;   ENDIF 
    250350      !      !--- dynamics without ridging/rafting and corr   (rheology + advection) 
    251351      IF( ln_dynRHGADV ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV    ;   ENDIF 
    252       !      !--- advection only with prescribed ice velocities (from namelist) 
    253       IF( ln_dynADV    ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV       ;   ENDIF 
     352      !      !--- advection 1D only - test case from Schar & Smolarkiewicz 1996 
     353      IF( ln_dynADV1D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV1D     ;   ENDIF 
     354      !      !--- advection 2D only with prescribed ice velocities (from namelist) 
     355      IF( ln_dynADV2D  ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV2D     ;   ENDIF 
    254356      ! 
    255357      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' ) 
     
    261363      ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  strong-slip' 
    262364      ENDIF 
    263       !                                      !--- NO Landfast ice : set to zero once for all 
    264       IF( .NOT.ln_landfast )   tau_icebfr(:,:) = 0._wp 
     365      !                                      !--- Landfast ice 
     366      IF( .NOT.ln_landfast_L16 .AND. .NOT.ln_landfast_home )   tau_icebfr(:,:) = 0._wp 
     367      ! 
     368      IF ( ln_landfast_L16 .AND. ln_landfast_home ) THEN 
     369         CALL ctl_stop( 'ice_dyn_init: choose one and only one landfast parameterization (ln_landfast_L16 or ln_landfast_home)' ) 
     370      ENDIF 
    265371      ! 
    266372      CALL ice_dyn_rdgrft_init          ! set ice ridging/rafting parameters 
Note: See TracChangeset for help on using the changeset viewer.