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 5053 for branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 – NEMO

Ignore:
Timestamp:
2015-02-03T18:11:02+01:00 (9 years ago)
Author:
clem
Message:

LIM3 cleaning continues. No change in the physics besides the introduction of the monocategory sea ice

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r5051 r5053  
    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 
    3535   !!---------------------------------------------------------------------- 
    3636   !!   'key_lim3'                                      LIM3 sea-ice model 
    37    !!---------------------------------------------------------------------- 
    38    !!   lim_var_agg       :  
    39    !!   lim_var_glo2eqv   : 
    40    !!   lim_var_eqv2glo   : 
    41    !!   lim_var_salprof   :  
    42    !!   lim_var_salprof1d : 
    43    !!   lim_var_bv        : 
    4437   !!---------------------------------------------------------------------- 
    4538   USE par_oce        ! ocean parameters 
     
    5851   PRIVATE 
    5952 
    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    ! 
     53   PUBLIC   lim_var_agg           
     54   PUBLIC   lim_var_glo2eqv       
     55   PUBLIC   lim_var_eqv2glo       
     56   PUBLIC   lim_var_salprof       
     57   PUBLIC   lim_var_icetm         
     58   PUBLIC   lim_var_bv            
     59   PUBLIC   lim_var_salprof1d     
    6760   PUBLIC   lim_var_zapsmall 
     61   PUBLIC   lim_var_itd 
    6862 
    6963   !!---------------------------------------------------------------------- 
    70    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     64   !! NEMO/LIM3 3.5 , UCL - NEMO Consortium (2011) 
    7165   !! $Id$ 
    7266   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    130124            DO jj = 1, jpj 
    131125               DO ji = 1, jpi 
    132                   et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                       ! snow heat content 
     126                  et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                           ! snow heat content 
    133127                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )  
    134128                  smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * rswitch   ! ice salinity 
     
    192186      ! Ice temperatures 
    193187      !------------------- 
    194 !CDIR NOVERRCHK 
    195       DO jl = 1, jpl 
    196 !CDIR NOVERRCHK 
     188      DO jl = 1, jpl 
    197189         DO jk = 1, nlay_i 
    198 !CDIR NOVERRCHK 
    199             DO jj = 1, jpj 
    200 !CDIR NOVERRCHK 
     190            DO jj = 1, jpj 
    201191               DO ji = 1, jpi 
    202192                  !                                                              ! Energy of melting q(S,T) [J.m-3] 
    203193                  rswitch   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) )     ! rswitch = 0 if no ice and 1 if yes 
    204194                  zq_i    = rswitch * e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp)  
    205                   zq_i    = zq_i * unit_fac                             !convert units 
    206                   ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt                       ! Ice layer melt temperature 
     195                  zq_i    = zq_i * unit_fac                             ! convert units 
     196                  ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt              ! Ice layer melt temperature 
    207197                  ! 
    208198                  zaaa       =  cpic                  ! Conversion q(S,T) -> T (second order equation) 
     
    289279      !!              alpha-weighted linear combination of s_inf and s_zero 
    290280      !! 
    291       !! ** References : Vancoppenolle et al., 2007 (in preparation) 
     281      !! ** References : Vancoppenolle et al., 2007 
    292282      !!------------------------------------------------------------------ 
    293283      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index 
     
    368358         ! 
    369359         DO jl = 1, jpl 
    370 !CDIR NOVERRCHK 
    371360            DO jk = 1, nlay_i 
    372361               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp) 
     
    451440      ! 
    452441      INTEGER  ::   ji, jk    ! dummy loop indices 
    453       INTEGER  ::   ii, ij  ! local integers 
     442      INTEGER  ::   ii, ij    ! local integers 
    454443      REAL(wp) ::   dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal   ! local scalars 
    455444      REAL(wp) ::   zalpha, zswi0, zswi01, zswibal, zs_zero              !   -      - 
     
    483472         dummy_fac2 = 1._wp / REAL(nlay_i,wp) 
    484473 
    485 !CDIR NOVERRCHK 
    486474         DO jk = 1, nlay_i 
    487 !CDIR NOVERRCHK 
    488475            DO ji = kideb, kiut 
    489476               ii =  MOD( npb(ji) - 1 , jpi ) + 1 
     
    502489               ! weighting the profile 
    503490               s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji) 
    504             END DO ! ji 
    505          END DO ! jk 
    506  
    507       ENDIF ! num_sal 
     491            END DO  
     492         END DO  
     493 
     494      ENDIF  
    508495 
    509496      !------------------------------------------------------- 
     
    515502         sm_i_1d(:) = 2.30_wp 
    516503         ! 
    517 !CDIR NOVERRCHK 
    518504         DO jk = 1, nlay_i 
    519505            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp) 
     
    603589      ! 
    604590   END SUBROUTINE lim_var_zapsmall 
     591 
     592   SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i ) 
     593      !!------------------------------------------------------------------ 
     594      !!                ***  ROUTINE lim_var_itd   *** 
     595      !! 
     596      !! ** Purpose :  converting 1-cat ice to multiple ice categories 
     597      !! 
     598      !!                  ice thickness distribution follows a gaussian law 
     599      !!               around the concentration of the most likely ice thickness 
     600      !!                           (similar as limistate.F90) 
     601      !! 
     602      !! ** Method:   Iterative procedure 
     603      !!                 
     604      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian 
     605      !! 
     606      !!               2) Check whether the distribution conserves area and volume, positivity and 
     607      !!                  category boundaries 
     608      !!               
     609      !!               3) If not (input ice is too thin), the last category is empty and 
     610      !!                  the number of categories is reduced (jpl-1) 
     611      !! 
     612      !!               4) Iterate until ok (SUM(itest(:) = 4) 
     613      !! 
     614      !! ** Arguments : zhti: 1-cat ice thickness 
     615      !!                zhts: 1-cat snow depth 
     616      !!                zai : 1-cat ice concentration 
     617      !! 
     618      !! ** Output    : jpl-cat  
     619      !! 
     620      !!  (Example of application: BDY forcings when input are cell averaged)   
     621      !! 
     622      !!------------------------------------------------------------------- 
     623      !! History : LIM3.5 - 2012    (M. Vancoppenolle)  Original code 
     624      !!                    2014    (C. Rousset)        Rewriting 
     625      !!------------------------------------------------------------------- 
     626      !! Local variables 
     627      INTEGER  :: ji, jk, jl             ! dummy loop indices 
     628      INTEGER  :: ijpij, i_fill, jl0   
     629      REAL(wp) :: zarg, zV, zconv, zdh 
     630      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables 
     631      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables 
     632      INTEGER , POINTER, DIMENSION(:)         ::   itest 
     633  
     634      CALL wrk_alloc( 4, itest ) 
     635      !-------------------------------------------------------------------- 
     636      ! initialisation of variables 
     637      !-------------------------------------------------------------------- 
     638      ijpij = SIZE(zhti,1) 
     639      zht_i(1:ijpij,1:jpl) = 0._wp 
     640      zht_s(1:ijpij,1:jpl) = 0._wp 
     641      za_i (1:ijpij,1:jpl) = 0._wp 
     642 
     643      ! ---------------------------------------- 
     644      ! distribution over the jpl ice categories 
     645      ! ---------------------------------------- 
     646      DO ji = 1, ijpij 
     647          
     648         IF( zhti(ji) > 0._wp ) THEN 
     649 
     650         ! initialisation of tests 
     651         itest(:)  = 0 
     652          
     653         i_fill = jpl + 1                                             !==================================== 
     654         DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories   
     655            ! iteration                                               !==================================== 
     656            i_fill = i_fill - 1 
     657             
     658            ! initialisation of ice variables for each try 
     659            zht_i(ji,1:jpl) = 0._wp 
     660            za_i (ji,1:jpl) = 0._wp 
     661             
     662            ! *** case very thin ice: fill only category 1 
     663            IF ( i_fill == 1 ) THEN 
     664               zht_i(ji,1) = zhti(ji) 
     665               za_i (ji,1) = zai (ji) 
     666 
     667            ! *** case ice is thicker: fill categories >1 
     668            ELSE 
     669 
     670               ! Fill ice thicknesses except the last one (i_fill) by hmean  
     671               DO jl = 1, i_fill - 1 
     672                  zht_i(ji,jl) = hi_mean(jl) 
     673               END DO 
     674                
     675               ! find which category (jl0) the input ice thickness falls into 
     676               jl0 = i_fill 
     677               DO jl = 1, i_fill 
     678                  IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 
     679                     jl0 = jl 
     680           CYCLE 
     681                  ENDIF 
     682               END DO 
     683                
     684               ! Concentrations in the (i_fill-1) categories  
     685               za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl)) 
     686               DO jl = 1, i_fill - 1 
     687                  IF ( jl == jl0 ) CYCLE 
     688                  zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 
     689                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2) 
     690               END DO 
     691                
     692               ! Concentration in the last (i_fill) category 
     693               za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) ) 
     694                
     695               ! Ice thickness in the last (i_fill) category 
     696               zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) ) 
     697               zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill)  
     698                
     699            ENDIF ! case ice is thick or thin 
     700             
     701            !--------------------- 
     702            ! Compatibility tests 
     703            !---------------------  
     704            ! Test 1: area conservation 
     705            zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) ) 
     706            IF ( zconv < epsi06 ) itest(1) = 1 
     707             
     708            ! Test 2: volume conservation 
     709            zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) ) 
     710            IF ( zconv < epsi06 ) itest(2) = 1 
     711             
     712            ! Test 3: thickness of the last category is in-bounds ? 
     713            IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 
     714             
     715            ! Test 4: positivity of ice concentrations 
     716            itest(4) = 1 
     717            DO jl = 1, i_fill 
     718               IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 
     719            END DO             
     720                                                           !============================ 
     721         END DO                                            ! end iteration on categories 
     722                                                           !============================ 
     723         ENDIF ! if zhti > 0 
     724      END DO ! i loop 
     725 
     726      ! ------------------------------------------------ 
     727      ! Adding Snow in each category where za_i is not 0 
     728      ! ------------------------------------------------  
     729      DO jl = 1, jpl 
     730         DO ji = 1, ijpij 
     731            IF( za_i(ji,jl) > 0._wp ) THEN 
     732               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) ) 
     733               ! In case snow load is in excess that would lead to transformation from snow to ice 
     734               ! Then, transfer the snow excess into the ice (different from limthd_dh) 
     735               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 )  
     736               ! recompute ht_i, ht_s avoiding out of bounds values 
     737               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh ) 
     738               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic / rhosn ) 
     739            ENDIF 
     740         ENDDO 
     741      ENDDO 
     742 
     743      CALL wrk_dealloc( 4, itest ) 
     744      ! 
     745    END SUBROUTINE lim_var_itd 
     746 
    605747 
    606748#else 
     
    623765   SUBROUTINE lim_var_zapsmall 
    624766   END SUBROUTINE lim_var_zapsmall 
     767   SUBROUTINE lim_var_itd 
     768   END SUBROUTINE lim_var_itd 
    625769#endif 
    626770 
Note: See TracChangeset for help on using the changeset viewer.