Ignore:
Timestamp:
2015-07-15T17:46:12+02:00 (5 years ago)
Author:
andrewryan
Message:

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r5034 r5600  
    3030   !!====================================================================== 
    3131   !! History :   -   ! 2006-01 (M. Vancoppenolle) Original code 
    32    !!            4.0  ! 2011-02 (G. Madec) dynamical allocation 
     32   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation 
    3333   !!---------------------------------------------------------------------- 
    3434#if defined key_lim3 
     
    3636   !!   'key_lim3'                                      LIM3 sea-ice model 
    3737   !!---------------------------------------------------------------------- 
    38    !!   lim_var_agg       :  
    39    !!   lim_var_glo2eqv   : 
    40    !!   lim_var_eqv2glo   : 
    41    !!   lim_var_salprof   :  
    42    !!   lim_var_salprof1d : 
    43    !!   lim_var_bv        : 
    44    !!---------------------------------------------------------------------- 
    4538   USE par_oce        ! ocean parameters 
    4639   USE phycst         ! physical constants (ocean directory)  
    4740   USE sbc_oce        ! Surface boundary condition: ocean fields 
    4841   USE ice            ! ice variables 
    49    USE par_ice        ! ice parameters 
    5042   USE thd_ice        ! ice variables (thermodynamics) 
    5143   USE dom_ice        ! ice domain 
     
    5850   PRIVATE 
    5951 
    60    PUBLIC   lim_var_agg          ! 
    61    PUBLIC   lim_var_glo2eqv      ! 
    62    PUBLIC   lim_var_eqv2glo      ! 
    63    PUBLIC   lim_var_salprof      ! 
    64    PUBLIC   lim_var_icetm        ! 
    65    PUBLIC   lim_var_bv           ! 
    66    PUBLIC   lim_var_salprof1d    ! 
     52   PUBLIC   lim_var_agg           
     53   PUBLIC   lim_var_glo2eqv       
     54   PUBLIC   lim_var_eqv2glo       
     55   PUBLIC   lim_var_salprof       
     56   PUBLIC   lim_var_icetm         
     57   PUBLIC   lim_var_bv            
     58   PUBLIC   lim_var_salprof1d     
     59   PUBLIC   lim_var_zapsmall 
     60   PUBLIC   lim_var_itd 
    6761 
    6862   !!---------------------------------------------------------------------- 
    69    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     63   !! NEMO/LIM3 3.5 , UCL - NEMO Consortium (2011) 
    7064   !! $Id$ 
    7165   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    129123            DO jj = 1, jpj 
    130124               DO ji = 1, jpi 
    131                   et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                       ! snow heat content 
    132                   rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )  
    133                   smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * rswitch   ! ice salinity 
    134                   rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )  
    135                   ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi10 ) * rswitch   ! ice age 
     125                  et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                           ! snow heat content 
     126                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi20 ) )  
     127                  smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi20 ) * rswitch   ! ice salinity 
     128                  rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi20 ) )  
     129                  ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi20 ) * rswitch   ! ice age 
    136130               END DO 
    137131            END DO 
     
    167161         DO jj = 1, jpj 
    168162            DO ji = 1, jpi 
    169                rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !0 if no ice and 1 if yes 
    170                ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 
    171                ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 
    172                o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 
    173             END DO 
    174          END DO 
    175       END DO 
    176  
    177       IF(  num_sal == 2  )THEN 
     163               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes 
     164               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
     165               ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
     166               o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
     167            END DO 
     168         END DO 
     169      END DO 
     170 
     171      IF(  nn_icesal == 2  )THEN 
    178172         DO jl = 1, jpl 
    179173            DO jj = 1, jpj 
    180174               DO ji = 1, jpi 
    181                   rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !0 if no ice and 1 if yes 
    182                   sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi10 ) * rswitch 
     175                  rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes 
     176                  sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * rswitch 
     177                  !                                      ! bounding salinity 
     178                  sm_i(ji,jj,jl) = MAX( sm_i(ji,jj,jl), rn_simin ) 
    183179               END DO 
    184180            END DO 
     
    191187      ! Ice temperatures 
    192188      !------------------- 
    193 !CDIR NOVERRCHK 
    194       DO jl = 1, jpl 
    195 !CDIR NOVERRCHK 
     189      DO jl = 1, jpl 
    196190         DO jk = 1, nlay_i 
    197 !CDIR NOVERRCHK 
    198191            DO jj = 1, jpj 
    199 !CDIR NOVERRCHK 
    200192               DO ji = 1, jpi 
    201193                  !                                                              ! Energy of melting q(S,T) [J.m-3] 
    202                   rswitch   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) )     ! rswitch = 0 if no ice and 1 if yes 
    203                   zq_i    = rswitch * e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp)  
    204                   zq_i    = zq_i * unit_fac                             !convert units 
    205                   ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt                       ! Ice layer melt temperature 
     194                  rswitch = MAX( 0.0 , SIGN( 1.0 , v_i(ji,jj,jl) - epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes 
     195                  zq_i    = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp)  
     196                  ztmelts = -tmut * s_i(ji,jj,jk,jl) + rt0              ! Ice layer melt temperature 
    206197                  ! 
    207198                  zaaa       =  cpic                  ! Conversion q(S,T) -> T (second order equation) 
    208                   zbbb       =  ( rcp - cpic ) * ( ztmelts - rtt ) + zq_i / rhoic - lfus 
    209                   zccc       =  lfus * (ztmelts-rtt) 
     199                  zbbb       =  ( rcp - cpic ) * ( ztmelts - rt0 ) + zq_i * r1_rhoic - lfus 
     200                  zccc       =  lfus * (ztmelts-rt0) 
    210201                  zdiscrim   =  SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) ) 
    211                   t_i(ji,jj,jk,jl) = rtt + rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa ) 
    212                   t_i(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15_wp, t_i(ji,jj,jk,jl) ) )       ! 100-rtt < t_i < rtt 
     202                  t_i(ji,jj,jk,jl) = rt0 + rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa ) 
     203                  t_i(ji,jj,jk,jl) = MIN( ztmelts, MAX( rt0 - 100._wp, t_i(ji,jj,jk,jl) ) )  ! -100 < t_i < ztmelts 
    213204               END DO 
    214205            END DO 
     
    226217               DO ji = 1, jpi 
    227218                  !Energy of melting q(S,T) [J.m-3] 
    228                   rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) )     ! rswitch = 0 if no ice and 1 if yes 
    229                   zq_s  = rswitch * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp) 
    230                   zq_s  = zq_s * unit_fac                                    ! convert units 
     219                  rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes 
     220                  zq_s  = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp) 
    231221                  ! 
    232                   t_s(ji,jj,jk,jl) = rtt + rswitch * ( - zfac1 * zq_s + zfac2 ) 
    233                   t_s(ji,jj,jk,jl) = MIN( rtt, MAX( 173.15, t_s(ji,jj,jk,jl) ) )     ! 100-rtt < t_i < rtt 
     222                  t_s(ji,jj,jk,jl) = rt0 + rswitch * ( - zfac1 * zq_s + zfac2 ) 
     223                  t_s(ji,jj,jk,jl) = MIN( rt0, MAX( rt0 - 100._wp , t_s(ji,jj,jk,jl) ) )     ! -100 < t_s < rt0 
    234224               END DO 
    235225            END DO 
     
    240230      ! Mean temperature 
    241231      !------------------- 
     232      vt_i (:,:) = 0._wp 
     233      DO jl = 1, jpl 
     234         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 
     235      END DO 
     236 
    242237      tm_i(:,:) = 0._wp 
    243238      DO jl = 1, jpl 
     
    245240            DO jj = 1, jpj 
    246241               DO ji = 1, jpi 
    247                   rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  ) 
    248                   tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   & 
    249                      &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 )  ) 
    250                END DO 
    251             END DO 
    252          END DO 
    253       END DO 
     242                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 
     243                  tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl)  & 
     244                     &            / MAX( vt_i(ji,jj) , epsi10 ) 
     245               END DO 
     246            END DO 
     247         END DO 
     248      END DO 
     249      tm_i = tm_i + rt0 
    254250      ! 
    255251   END SUBROUTINE lim_var_glo2eqv 
     
    270266      v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:) 
    271267      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 
    272       oa_i (:,:,:) = o_i (:,:,:) * a_i(:,:,:) 
    273268      ! 
    274269   END SUBROUTINE lim_var_eqv2glo 
     
    281276      !! ** Purpose :   computes salinity profile in function of bulk salinity      
    282277      !! 
    283       !! ** Method  : If bulk salinity greater than s_i_1,  
     278      !! ** Method  : If bulk salinity greater than zsi1,  
    284279      !!              the profile is assumed to be constant (S_inf) 
    285       !!              If bulk salinity lower than s_i_0, 
     280      !!              If bulk salinity lower than zsi0, 
    286281      !!              the profile is linear with 0 at the surface (S_zero) 
    287       !!              If it is between s_i_0 and s_i_1, it is a 
     282      !!              If it is between zsi0 and zsi1, it is a 
    288283      !!              alpha-weighted linear combination of s_inf and s_zero 
    289284      !! 
    290       !! ** References : Vancoppenolle et al., 2007 (in preparation) 
     285      !! ** References : Vancoppenolle et al., 2007 
    291286      !!------------------------------------------------------------------ 
    292287      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index 
    293       REAL(wp) ::   dummy_fac0, dummy_fac1, dummy_fac, zsal      ! local scalar 
    294       REAL(wp) ::   zswi0, zswi01, zswibal, zargtemp , zs_zero   !   -      - 
    295       REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_slope_s, zalpha   ! 3D pointer 
     288      REAL(wp) ::   zfac0, zfac1, zsal 
     289      REAL(wp) ::   zswi0, zswi01, zargtemp , zs_zero    
     290      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_slope_s, zalpha 
     291      REAL(wp), PARAMETER :: zsi0 = 3.5_wp 
     292      REAL(wp), PARAMETER :: zsi1 = 4.5_wp 
    296293      !!------------------------------------------------------------------ 
    297294 
     
    301298      ! Vertically constant, constant in time 
    302299      !--------------------------------------- 
    303       IF(  num_sal == 1  )   s_i(:,:,:,:) = bulk_sal 
     300      IF(  nn_icesal == 1  )   s_i(:,:,:,:) = rn_icesal 
    304301 
    305302      !----------------------------------- 
    306303      ! Salinity profile, varying in time 
    307304      !----------------------------------- 
    308       IF(  num_sal == 2  ) THEN 
     305      IF(  nn_icesal == 2  ) THEN 
    309306         ! 
    310307         DO jk = 1, nlay_i 
     
    315312            DO jj = 1, jpj 
    316313               DO ji = 1, jpi 
    317                   z_slope_s(ji,jj,jl) = 2._wp * sm_i(ji,jj,jl) / MAX( epsi10 , ht_i(ji,jj,jl) ) 
    318                END DO 
    319             END DO 
    320          END DO 
    321          ! 
    322          dummy_fac0 = 1._wp / ( s_i_0 - s_i_1 )       ! Weighting factor between zs_zero and zs_inf 
    323          dummy_fac1 = s_i_1 / ( s_i_1 - s_i_0 ) 
     314                  rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i(ji,jj,jl) - epsi20 ) ) 
     315                  z_slope_s(ji,jj,jl) = rswitch * 2._wp * sm_i(ji,jj,jl) / MAX( epsi20 , ht_i(ji,jj,jl) ) 
     316               END DO 
     317            END DO 
     318         END DO 
     319         ! 
     320         zfac0 = 1._wp / ( zsi0 - zsi1 )       ! Weighting factor between zs_zero and zs_inf 
     321         zfac1 = zsi1  / ( zsi1 - zsi0 ) 
    324322         ! 
    325323         zalpha(:,:,:) = 0._wp 
     
    327325            DO jj = 1, jpj 
    328326               DO ji = 1, jpi 
    329                   ! zswi0 = 1 if sm_i le s_i_0 and 0 otherwise 
    330                   zswi0  = MAX( 0._wp   , SIGN( 1._wp  , s_i_0 - sm_i(ji,jj,jl) ) )  
    331                   ! zswi01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws  
    332                   zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp   , SIGN( 1._wp  , s_i_1 - sm_i(ji,jj,jl) ) )  
    333                   ! If 2.sm_i GE sss_m then zswibal = 1 
     327                  ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise 
     328                  zswi0  = MAX( 0._wp   , SIGN( 1._wp  , zsi0 - sm_i(ji,jj,jl) ) )  
     329                  ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws  
     330                  zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp   , SIGN( 1._wp  , zsi1 - sm_i(ji,jj,jl) ) )  
     331                  ! If 2.sm_i GE sss_m then rswitch = 1 
    334332                  ! this is to force a constant salinity profile in the Baltic Sea 
    335                   zswibal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) ) 
    336                   zalpha(ji,jj,jl) = zswi0  + zswi01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 ) 
    337                   zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - zswibal ) 
    338                END DO 
    339             END DO 
    340          END DO 
    341  
    342          dummy_fac = 1._wp / REAL( nlay_i )                   ! Computation of the profile 
     333                  rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) ) 
     334                  zalpha(ji,jj,jl) = zswi0  + zswi01 * ( sm_i(ji,jj,jl) * zfac0 + zfac1 ) 
     335                  zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - rswitch ) 
     336               END DO 
     337            END DO 
     338         END DO 
     339 
     340         ! Computation of the profile 
    343341         DO jl = 1, jpl 
    344342            DO jk = 1, nlay_i 
     
    346344                  DO ji = 1, jpi 
    347345                     !                                      ! linear profile with 0 at the surface 
    348                      zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * dummy_fac 
     346                     zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i 
    349347                     !                                      ! weighting the profile 
    350348                     s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl) 
    351                   END DO ! ji 
    352                END DO ! jj 
    353             END DO ! jk 
    354          END DO ! jl 
    355          ! 
    356       ENDIF ! num_sal 
     349                     !                                      ! bounding salinity 
     350                     s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( s_i(ji,jj,jk,jl), rn_simin ) ) 
     351                  END DO 
     352               END DO 
     353            END DO 
     354         END DO 
     355         ! 
     356      ENDIF ! nn_icesal 
    357357 
    358358      !------------------------------------------------------- 
     
    360360      !------------------------------------------------------- 
    361361 
    362       IF(  num_sal == 3  ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 
     362      IF(  nn_icesal == 3  ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 
    363363         ! 
    364364         sm_i(:,:,:) = 2.30_wp 
    365365         ! 
    366366         DO jl = 1, jpl 
    367 !CDIR NOVERRCHK 
    368367            DO jk = 1, nlay_i 
    369                zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp) 
     368               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i 
    370369               zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  ) 
    371370               s_i(:,:,jk,jl) =  zsal 
     
    373372         END DO 
    374373         ! 
    375       ENDIF ! num_sal 
     374      ENDIF ! nn_icesal 
    376375      ! 
    377376      CALL wrk_dealloc( jpi, jpj, jpl, z_slope_s, zalpha ) 
     
    390389 
    391390      ! Mean sea ice temperature 
     391      vt_i (:,:) = 0._wp 
     392      DO jl = 1, jpl 
     393         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 
     394      END DO 
     395 
    392396      tm_i(:,:) = 0._wp 
    393397      DO jl = 1, jpl 
     
    395399            DO jj = 1, jpj 
    396400               DO ji = 1, jpi 
    397                   rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  ) 
    398                   tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   & 
    399                      &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 )  ) 
    400                END DO 
    401             END DO 
    402          END DO 
    403       END DO 
     401                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 
     402                  tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl)  & 
     403                     &            / MAX( vt_i(ji,jj) , epsi10 ) 
     404               END DO 
     405            END DO 
     406         END DO 
     407      END DO 
     408      tm_i = tm_i + rt0 
    404409 
    405410   END SUBROUTINE lim_var_icetm 
     
    420425      !!------------------------------------------------------------------ 
    421426      ! 
     427      vt_i (:,:) = 0._wp 
     428      DO jl = 1, jpl 
     429         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 
     430      END DO 
     431 
    422432      bv_i(:,:) = 0._wp 
    423433      DO jl = 1, jpl 
     
    425435            DO jj = 1, jpj 
    426436               DO ji = 1, jpi 
    427                   rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) )  ) 
    428                   zbvi  = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 )   & 
    429                      &                   * v_i(ji,jj,jl)    / REAL(nlay_i,wp) 
    430                   rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  ) 
    431                   bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi  / MAX( vt_i(ji,jj) , epsi10 ) 
     437                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) )  ) 
     438                  zbvi  = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )   & 
     439                     &                   * v_i(ji,jj,jl) * r1_nlay_i 
     440                  rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi20 ) )  ) 
     441                  bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi  / MAX( vt_i(ji,jj) , epsi20 ) 
    432442               END DO 
    433443            END DO 
     
    448458      ! 
    449459      INTEGER  ::   ji, jk    ! dummy loop indices 
    450       INTEGER  ::   ii, ij  ! local integers 
    451       REAL(wp) ::   dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal   ! local scalars 
    452       REAL(wp) ::   zalpha, zswi0, zswi01, zswibal, zs_zero              !   -      - 
     460      INTEGER  ::   ii, ij    ! local integers 
     461      REAL(wp) ::   zfac0, zfac1, zargtemp, zsal   ! local scalars 
     462      REAL(wp) ::   zalpha, zswi0, zswi01, zs_zero              !   -      - 
    453463      ! 
    454464      REAL(wp), POINTER, DIMENSION(:) ::   z_slope_s 
     465      REAL(wp), PARAMETER :: zsi0 = 3.5_wp 
     466      REAL(wp), PARAMETER :: zsi1 = 4.5_wp 
    455467      !!--------------------------------------------------------------------- 
    456468 
     
    460472      ! Vertically constant, constant in time 
    461473      !--------------------------------------- 
    462       IF( num_sal == 1 )   s_i_1d(:,:) = bulk_sal 
     474      IF( nn_icesal == 1 )   s_i_1d(:,:) = rn_icesal 
    463475 
    464476      !------------------------------------------------------ 
     
    466478      !------------------------------------------------------ 
    467479 
    468       IF(  num_sal == 2  ) THEN 
     480      IF(  nn_icesal == 2  ) THEN 
    469481         ! 
    470482         DO ji = kideb, kiut          ! Slope of the linear profile zs_zero 
    471             z_slope_s(ji) = 2._wp * sm_i_1d(ji) / MAX( epsi10 , ht_i_1d(ji) ) 
     483            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) ) 
     484            z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) ) 
    472485         END DO 
    473486 
    474487         ! Weighting factor between zs_zero and zs_inf 
    475488         !--------------------------------------------- 
    476          dummy_fac0 = 1._wp / ( s_i_0 - s_i_1 ) 
    477          dummy_fac1 = s_i_1 / ( s_i_1 - s_i_0 ) 
    478          dummy_fac2 = 1._wp / REAL(nlay_i,wp) 
    479  
    480 !CDIR NOVERRCHK 
     489         zfac0 = 1._wp / ( zsi0 - zsi1 ) 
     490         zfac1 = zsi1 / ( zsi1 - zsi0 ) 
    481491         DO jk = 1, nlay_i 
    482 !CDIR NOVERRCHK 
    483492            DO ji = kideb, kiut 
    484493               ii =  MOD( npb(ji) - 1 , jpi ) + 1 
    485494               ij =     ( npb(ji) - 1 ) / jpi + 1 
    486                ! zswi0 = 1 if sm_i le s_i_0 and 0 otherwise 
    487                zswi0  = MAX( 0._wp , SIGN( 1._wp  , s_i_0 - sm_i_1d(ji) ) )  
    488                ! zswi01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws  
    489                zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_1d(ji) ) )  
    490                ! if 2.sm_i GE sss_m then zswibal = 1 
     495               ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise 
     496               zswi0  = MAX( 0._wp , SIGN( 1._wp  , zsi0 - sm_i_1d(ji) ) )  
     497               ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws  
     498               zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) )  
     499               ! if 2.sm_i GE sss_m then rswitch = 1 
    491500               ! this is to force a constant salinity profile in the Baltic Sea 
    492                zswibal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) ) 
     501               rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) ) 
    493502               ! 
    494                zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * dummy_fac0 + dummy_fac1 )  ) * ( 1.0 - zswibal ) 
     503               zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 )  ) * ( 1._wp - rswitch ) 
    495504               ! 
    496                zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * dummy_fac2 
     505               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i 
    497506               ! weighting the profile 
    498507               s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji) 
    499             END DO ! ji 
    500          END DO ! jk 
    501  
    502       ENDIF ! num_sal 
     508               ! bounding salinity 
     509               s_i_1d(ji,jk) = MIN( rn_simax, MAX( s_i_1d(ji,jk), rn_simin ) ) 
     510            END DO  
     511         END DO  
     512 
     513      ENDIF  
    503514 
    504515      !------------------------------------------------------- 
     
    506517      !------------------------------------------------------- 
    507518 
    508       IF( num_sal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 
     519      IF( nn_icesal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 
    509520         ! 
    510521         sm_i_1d(:) = 2.30_wp 
    511522         ! 
    512 !CDIR NOVERRCHK 
    513523         DO jk = 1, nlay_i 
    514             zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp) 
    515             zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) ) 
     524            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i 
     525            zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) ) 
    516526            DO ji = kideb, kiut 
    517527               s_i_1d(ji,jk) = zsal 
     
    524534      ! 
    525535   END SUBROUTINE lim_var_salprof1d 
     536 
     537   SUBROUTINE lim_var_zapsmall 
     538      !!------------------------------------------------------------------- 
     539      !!                   ***  ROUTINE lim_var_zapsmall *** 
     540      !! 
     541      !! ** Purpose :   Remove too small sea ice areas and correct fluxes 
     542      !! 
     543      !! history : LIM3.5 - 01-2014 (C. Rousset) original code 
     544      !!------------------------------------------------------------------- 
     545      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices 
     546      REAL(wp) ::   zsal, zvi, zvs, zei, zes 
     547      !!------------------------------------------------------------------- 
     548      at_i (:,:) = 0._wp 
     549      DO jl = 1, jpl 
     550         at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
     551      END DO 
     552 
     553      DO jl = 1, jpl 
     554 
     555         !----------------------------------------------------------------- 
     556         ! Zap ice energy and use ocean heat to melt ice 
     557         !----------------------------------------------------------------- 
     558         DO jk = 1, nlay_i 
     559            DO jj = 1 , jpj 
     560               DO ji = 1 , jpi 
     561                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 
     562                  rswitch          = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj  ) - epsi10 ) ) * rswitch 
     563                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch 
     564                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  & 
     565                     &                                       / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch 
     566                  zei              = e_i(ji,jj,jk,jl) 
     567                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch 
     568                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch ) 
     569                  ! update exchanges with ocean 
     570                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0 
     571               END DO 
     572            END DO 
     573         END DO 
     574 
     575         DO jj = 1 , jpj 
     576            DO ji = 1 , jpi 
     577               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 
     578               rswitch = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj  ) - epsi10 ) ) * rswitch 
     579               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch 
     580               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  & 
     581                  &                              / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch 
     582               zsal = smv_i(ji,jj,  jl) 
     583               zvi  = v_i  (ji,jj,  jl) 
     584               zvs  = v_s  (ji,jj,  jl) 
     585               zes  = e_s  (ji,jj,1,jl) 
     586               !----------------------------------------------------------------- 
     587               ! Zap snow energy  
     588               !----------------------------------------------------------------- 
     589               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch ) 
     590               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch 
     591 
     592               !----------------------------------------------------------------- 
     593               ! zap ice and snow volume, add water and salt to ocean 
     594               !----------------------------------------------------------------- 
     595               ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj) 
     596               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * rswitch 
     597               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * rswitch 
     598               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * rswitch 
     599               t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch ) 
     600               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch 
     601               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch 
     602 
     603               ! update exchanges with ocean 
     604               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 
     605               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice 
     606               wfx_snw(ji,jj)  = wfx_snw(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice 
     607               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 
     608            END DO 
     609         END DO 
     610      END DO  
     611 
     612      ! to be sure that at_i is the sum of a_i(jl) 
     613      at_i (:,:) = 0._wp 
     614      DO jl = 1, jpl 
     615         at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
     616      END DO 
     617 
     618      ! open water = 1 if at_i=0 
     619      DO jj = 1, jpj 
     620         DO ji = 1, jpi 
     621            rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 
     622            ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj) 
     623         END DO 
     624      END DO 
     625 
     626      ! 
     627   END SUBROUTINE lim_var_zapsmall 
     628 
     629   SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i ) 
     630      !!------------------------------------------------------------------ 
     631      !!                ***  ROUTINE lim_var_itd   *** 
     632      !! 
     633      !! ** Purpose :  converting 1-cat ice to multiple ice categories 
     634      !! 
     635      !!                  ice thickness distribution follows a gaussian law 
     636      !!               around the concentration of the most likely ice thickness 
     637      !!                           (similar as limistate.F90) 
     638      !! 
     639      !! ** Method:   Iterative procedure 
     640      !!                 
     641      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian 
     642      !! 
     643      !!               2) Check whether the distribution conserves area and volume, positivity and 
     644      !!                  category boundaries 
     645      !!               
     646      !!               3) If not (input ice is too thin), the last category is empty and 
     647      !!                  the number of categories is reduced (jpl-1) 
     648      !! 
     649      !!               4) Iterate until ok (SUM(itest(:) = 4) 
     650      !! 
     651      !! ** Arguments : zhti: 1-cat ice thickness 
     652      !!                zhts: 1-cat snow depth 
     653      !!                zai : 1-cat ice concentration 
     654      !! 
     655      !! ** Output    : jpl-cat  
     656      !! 
     657      !!  (Example of application: BDY forcings when input are cell averaged)   
     658      !! 
     659      !!------------------------------------------------------------------- 
     660      !! History : LIM3.5 - 2012    (M. Vancoppenolle)  Original code 
     661      !!                    2014    (C. Rousset)        Rewriting 
     662      !!------------------------------------------------------------------- 
     663      !! Local variables 
     664      INTEGER  :: ji, jk, jl             ! dummy loop indices 
     665      INTEGER  :: ijpij, i_fill, jl0   
     666      REAL(wp) :: zarg, zV, zconv, zdh 
     667      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables 
     668      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables 
     669      INTEGER , POINTER, DIMENSION(:)         ::   itest 
     670  
     671      CALL wrk_alloc( 4, itest ) 
     672      !-------------------------------------------------------------------- 
     673      ! initialisation of variables 
     674      !-------------------------------------------------------------------- 
     675      ijpij = SIZE(zhti,1) 
     676      zht_i(1:ijpij,1:jpl) = 0._wp 
     677      zht_s(1:ijpij,1:jpl) = 0._wp 
     678      za_i (1:ijpij,1:jpl) = 0._wp 
     679 
     680      ! ---------------------------------------- 
     681      ! distribution over the jpl ice categories 
     682      ! ---------------------------------------- 
     683      DO ji = 1, ijpij 
     684          
     685         IF( zhti(ji) > 0._wp ) THEN 
     686 
     687         ! initialisation of tests 
     688         itest(:)  = 0 
     689          
     690         i_fill = jpl + 1                                             !==================================== 
     691         DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories   
     692            ! iteration                                               !==================================== 
     693            i_fill = i_fill - 1 
     694             
     695            ! initialisation of ice variables for each try 
     696            zht_i(ji,1:jpl) = 0._wp 
     697            za_i (ji,1:jpl) = 0._wp 
     698             
     699            ! *** case very thin ice: fill only category 1 
     700            IF ( i_fill == 1 ) THEN 
     701               zht_i(ji,1) = zhti(ji) 
     702               za_i (ji,1) = zai (ji) 
     703 
     704            ! *** case ice is thicker: fill categories >1 
     705            ELSE 
     706 
     707               ! Fill ice thicknesses except the last one (i_fill) by hmean  
     708               DO jl = 1, i_fill - 1 
     709                  zht_i(ji,jl) = hi_mean(jl) 
     710               END DO 
     711                
     712               ! find which category (jl0) the input ice thickness falls into 
     713               jl0 = i_fill 
     714               DO jl = 1, i_fill 
     715                  IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 
     716                     jl0 = jl 
     717           CYCLE 
     718                  ENDIF 
     719               END DO 
     720                
     721               ! Concentrations in the (i_fill-1) categories  
     722               za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl)) 
     723               DO jl = 1, i_fill - 1 
     724                  IF ( jl == jl0 ) CYCLE 
     725                  zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 
     726                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2) 
     727               END DO 
     728                
     729               ! Concentration in the last (i_fill) category 
     730               za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) ) 
     731                
     732               ! Ice thickness in the last (i_fill) category 
     733               zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) ) 
     734               zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill)  
     735                
     736            ENDIF ! case ice is thick or thin 
     737             
     738            !--------------------- 
     739            ! Compatibility tests 
     740            !---------------------  
     741            ! Test 1: area conservation 
     742            zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) ) 
     743            IF ( zconv < epsi06 ) itest(1) = 1 
     744             
     745            ! Test 2: volume conservation 
     746            zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) ) 
     747            IF ( zconv < epsi06 ) itest(2) = 1 
     748             
     749            ! Test 3: thickness of the last category is in-bounds ? 
     750            IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 
     751             
     752            ! Test 4: positivity of ice concentrations 
     753            itest(4) = 1 
     754            DO jl = 1, i_fill 
     755               IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 
     756            END DO             
     757                                                           !============================ 
     758         END DO                                            ! end iteration on categories 
     759                                                           !============================ 
     760         ENDIF ! if zhti > 0 
     761      END DO ! i loop 
     762 
     763      ! ------------------------------------------------ 
     764      ! Adding Snow in each category where za_i is not 0 
     765      ! ------------------------------------------------  
     766      DO jl = 1, jpl 
     767         DO ji = 1, ijpij 
     768            IF( za_i(ji,jl) > 0._wp ) THEN 
     769               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) ) 
     770               ! In case snow load is in excess that would lead to transformation from snow to ice 
     771               ! Then, transfer the snow excess into the ice (different from limthd_dh) 
     772               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 )  
     773               ! recompute ht_i, ht_s avoiding out of bounds values 
     774               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh ) 
     775               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn ) 
     776            ENDIF 
     777         ENDDO 
     778      ENDDO 
     779 
     780      CALL wrk_dealloc( 4, itest ) 
     781      ! 
     782    END SUBROUTINE lim_var_itd 
     783 
    526784 
    527785#else 
     
    542800   SUBROUTINE lim_var_salprof1d    ! Emtpy routines 
    543801   END SUBROUTINE lim_var_salprof1d 
     802   SUBROUTINE lim_var_zapsmall 
     803   END SUBROUTINE lim_var_zapsmall 
     804   SUBROUTINE lim_var_itd 
     805   END SUBROUTINE lim_var_itd 
    544806#endif 
    545807 
Note: See TracChangeset for help on using the changeset viewer.