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 4161 for branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 – NEMO

Ignore:
Timestamp:
2013-11-07T11:01:27+01:00 (10 years ago)
Author:
cetlod
Message:

dev_LOCEAN_2013 : merge in the 3rd dev branch dev_r4028_CNRS_LIM3, see ticket #1169

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r3625 r4161  
    6262   PUBLIC   lim_var_eqv2glo      ! 
    6363   PUBLIC   lim_var_salprof      ! 
     64   PUBLIC   lim_var_icetm        ! 
    6465   PUBLIC   lim_var_bv           ! 
    6566   PUBLIC   lim_var_salprof1d    ! 
    6667 
    67    REAL(wp) ::   eps20 = 1.e-20_wp   ! module constants 
    68    REAL(wp) ::   eps16 = 1.e-16_wp   !    -       - 
    69    REAL(wp) ::   eps13 = 1.e-13_wp   !    -       - 
    70    REAL(wp) ::   eps10 = 1.e-10_wp   !    -       - 
    71    REAL(wp) ::   eps06 = 1.e-06_wp   !    -       - 
     68   REAL(wp) ::   epsi20 = 1.e-20_wp   ! module constants 
     69   REAL(wp) ::   epsi16 = 1.e-16_wp   !    -       - 
     70   REAL(wp) ::   epsi13 = 1.e-13_wp   !    -       - 
     71   REAL(wp) ::   epsi10 = 1.e-10_wp   !    -       - 
     72   REAL(wp) ::   epsi06 = 1.e-06_wp   !    -       - 
    7273   REAL(wp) ::   zzero = 0.e0        !    -       - 
    7374   REAL(wp) ::   zone  = 1.e0        !    -       - 
    7475 
    7576   !!---------------------------------------------------------------------- 
    76    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     77   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    7778   !! $Id$ 
    7879   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    9798      ! 
    9899      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    99       REAL(wp) ::   zinda 
     100      REAL(wp) ::   zinda, zindb 
    100101      !!------------------------------------------------------------------ 
    101102 
     
    116117               at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 
    117118               ! 
    118                zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - 0.10 ) )  
    119                icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , eps16 ) * zinda  ! ice thickness 
     119               zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi16 ) )  
     120               icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi16 ) * zinda  ! ice thickness 
    120121            END DO 
    121122         END DO 
     
    137138            DO jj = 1, jpj 
    138139               DO ji = 1, jpi 
     140                  zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - epsi16 ) )  
     141                  zindb = MAX( zzero , SIGN( zone , at_i(ji,jj) - epsi16 ) )  
    139142                  et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                       ! snow heat content 
    140                   zinda = MAX( zzero , SIGN( zone , vt_i(ji,jj) - 0.10 ) )  
    141                   smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , eps13 ) * zinda   ! ice salinity 
    142                   zinda = MAX( zzero , SIGN( zone , at_i(ji,jj) - 0.10 ) )  
    143                   ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , eps13 ) * zinda   ! ice age 
     143                  smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi16 ) * zinda   ! ice salinity 
     144                  ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi16 ) * zindb   ! ice age 
    144145               END DO 
    145146            END DO 
     
    175176         DO jj = 1, jpj 
    176177            DO ji = 1, jpi 
    177                zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) ) )   !0 if no ice and 1 if yes 
    178                ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , eps10 ) * zindb 
    179                ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , eps10 ) * zindb 
    180                o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , eps10 ) * zindb 
     178               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !0 if no ice and 1 if yes 
     179               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb 
     180               ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb 
     181               o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb 
     182               a_i(ji,jj,jl) = a_i (ji,jj,jl) * zindb ! clem correction 
    181183            END DO 
    182184         END DO 
     
    187189            DO jj = 1, jpj 
    188190               DO ji = 1, jpi 
    189                   zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) ) )   !0 if no ice and 1 if yes 
    190                   sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , eps10 ) * zindb 
     191                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp,- a_i(ji,jj,jl) + epsi10 ) )   !0 if no ice and 1 if yes 
     192                  sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi10 ) * zindb 
    191193               END DO 
    192194            END DO 
     
    208210               DO ji = 1, jpi 
    209211                  !                                                              ! Energy of melting q(S,T) [J.m-3] 
    210                   zq_i    = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , eps06 ) * REAL(nlay_i,wp)  
     212                  zq_i    = e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi06 ) * REAL(nlay_i,wp)  
    211213                  zindb   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) ) )     ! zindb = 0 if no ice and 1 if yes 
    212214                  zq_i    = zq_i * unit_fac * zindb                              !convert units 
     
    234236               DO ji = 1, jpi 
    235237                  !Energy of melting q(S,T) [J.m-3] 
    236                   zq_s  = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , eps06 ) ) * REAL(nlay_s,wp) 
     238                  zq_s  = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * REAL(nlay_s,wp) 
    237239                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) ) )     ! zindb = 0 if no ice and 1 if yes 
    238240                  zq_s  = zq_s * unit_fac * zindb                                    ! convert units 
     
    253255            DO jj = 1, jpj 
    254256               DO ji = 1, jpi 
    255                   zindb = (  1._wp - MAX( 0._wp , SIGN( 1._wp , -a_i(ji,jj,jl) ) )  )   & 
    256                      &  * (  1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) ) )  ) 
    257                   tm_i(ji,jj) = tm_i(ji,jj) + t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   & 
    258                      &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , eps10 )  ) 
     257                  zindb = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  ) 
     258                  tm_i(ji,jj) = tm_i(ji,jj) + zindb * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   & 
     259                     &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 )  ) 
    259260               END DO 
    260261            END DO 
     
    337338               DO ji = 1, jpi 
    338339                  ! zind0 = 1 if sm_i le s_i_0 and 0 otherwise 
    339                   zind0  = MAX( 0.0   , SIGN( 1.0  , s_i_0 - sm_i(ji,jj,jl) ) )  
     340                  zind0  = MAX( 0._wp   , SIGN( 1._wp  , s_i_0 - sm_i(ji,jj,jl) ) )  
    340341                  ! zind01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws  
    341                   zind01 = ( 1.0 - zind0 ) * MAX( 0.0   , SIGN( 1.0  , s_i_1 - sm_i(ji,jj,jl) ) )  
     342                  zind01 = ( 1._wp - zind0 ) * MAX( 0._wp   , SIGN( 1._wp  , s_i_1 - sm_i(ji,jj,jl) ) )  
    342343                  ! If 2.sm_i GE sss_m then zindbal = 1 
    343                   zindbal = MAX( 0.0 , SIGN( 1.0 , 2. * sm_i(ji,jj,jl) - sss_m(ji,jj) ) ) 
    344                   zalpha(ji,jj,jl) = zind0  * 1.0 + zind01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 ) 
    345                   zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1.0 - zindbal ) 
    346                END DO 
    347             END DO 
    348          END DO 
    349          ! 
    350          dummy_fac = 1._wp / nlay_i                   ! Computation of the profile 
     344                  zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) ) 
     345                  zalpha(ji,jj,jl) = zind0  + zind01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 ) 
     346                  zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - zindbal ) 
     347               END DO 
     348            END DO 
     349         END DO 
     350 
     351         dummy_fac = 1._wp / REAL( nlay_i )                   ! Computation of the profile 
    351352         DO jl = 1, jpl 
    352353            DO jk = 1, nlay_i 
     
    388389 
    389390 
     391   SUBROUTINE lim_var_icetm 
     392      !!------------------------------------------------------------------ 
     393      !!                ***  ROUTINE lim_var_icetm *** 
     394      !! 
     395      !! ** Purpose :   computes mean sea ice temperature 
     396      !!------------------------------------------------------------------ 
     397      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
     398      REAL(wp) ::   zindb   !   -      - 
     399      !!------------------------------------------------------------------ 
     400 
     401      ! Mean sea ice temperature 
     402      tm_i(:,:) = 0._wp 
     403      DO jl = 1, jpl 
     404         DO jk = 1, nlay_i 
     405            DO jj = 1, jpj 
     406               DO ji = 1, jpi 
     407                  zindb = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  ) 
     408                  tm_i(ji,jj) = tm_i(ji,jj) + zindb * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   & 
     409                     &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 )  ) 
     410               END DO 
     411            END DO 
     412         END DO 
     413      END DO 
     414 
     415   END SUBROUTINE lim_var_icetm 
     416 
     417 
    390418   SUBROUTINE lim_var_bv 
    391419      !!------------------------------------------------------------------ 
     
    399427      !!------------------------------------------------------------------ 
    400428      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    401       REAL(wp) ::   zbvi, zindb      ! local scalars 
     429      REAL(wp) ::   zbvi, zinda, zindb      ! local scalars 
    402430      !!------------------------------------------------------------------ 
    403431      ! 
     
    407435            DO jj = 1, jpj 
    408436               DO ji = 1, jpi 
    409                   zindb = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl))) !0 if no ice and 1 if yes 
    410                   zbvi  = - zindb * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - 273.15 , eps13 )   & 
     437                  zinda = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi16 ) )  ) 
     438                  zindb = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi16 ) )  ) 
     439                  zbvi  = - zinda * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi16 )   & 
    411440                     &                   * v_i(ji,jj,jl)    / REAL(nlay_i,wp) 
    412                   bv_i(ji,jj) = bv_i(ji,jj) + zbvi  / MAX( vt_i(ji,jj) , eps13 ) 
     441                  bv_i(ji,jj) = bv_i(ji,jj) + zindb * zbvi  / MAX( vt_i(ji,jj) , epsi16 ) 
    413442               END DO 
    414443            END DO 
     
    429458      ! 
    430459      INTEGER  ::   ji, jk    ! dummy loop indices 
    431       INTEGER  ::   zji, zjj  ! local integers 
     460      INTEGER  ::   ii, ij  ! local integers 
    432461      REAL(wp) ::   dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal   ! local scalars 
    433462      REAL(wp) ::   zalpha, zind0, zind01, zindbal, zs_zero              !   -      - 
     
    463492!CDIR NOVERRCHK 
    464493            DO ji = kideb, kiut 
    465                zji =  MOD( npb(ji) - 1 , jpi ) + 1 
    466                zjj =     ( npb(ji) - 1 ) / jpi + 1 
     494               ii =  MOD( npb(ji) - 1 , jpi ) + 1 
     495               ij =     ( npb(ji) - 1 ) / jpi + 1 
    467496               ! zind0 = 1 if sm_i le s_i_0 and 0 otherwise 
    468497               zind0  = MAX( 0._wp , SIGN( 1._wp  , s_i_0 - sm_i_b(ji) ) )  
     
    470499               zind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_b(ji) ) )  
    471500               ! if 2.sm_i GE sss_m then zindbal = 1 
    472                zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_b(ji) - sss_m(zji,zjj) ) ) 
     501               zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_b(ji) - sss_m(ii,ij) ) ) 
    473502               ! 
    474503               zalpha = (  zind0 + zind01 * ( sm_i_b(ji) * dummy_fac0 + dummy_fac1 )  ) * ( 1.0 - zindbal ) 
Note: See TracChangeset for help on using the changeset viewer.