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 4332 for branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 – NEMO

Ignore:
Timestamp:
2013-12-11T15:38:42+01:00 (10 years ago)
Author:
clem
Message:

update LIM3 to fix remaining bugs. Now working in global and regional config.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r4220 r4332  
    3131   PUBLIC   lim_thd_dif   ! called by lim_thd 
    3232 
    33    REAL(wp) ::   epsi20 = 1e-20     ! constant values 
    34    REAL(wp) ::   epsi13 = 1e-13     ! constant values 
    35  
     33   REAL(wp) ::   epsi10      =  1.e-10_wp    ! 
    3634   !!---------------------------------------------------------------------- 
    3735   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     
    107105      INTEGER, DIMENSION(kiut) ::   numeqmax   ! reference number of bottom equation 
    108106      INTEGER, DIMENSION(kiut) ::   isnow      ! switch for presence (1) or absence (0) of snow 
    109       REAL(wp) ::   zeps      =  1.e-10_wp    ! 
    110107      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system 
    111108      REAL(wp) ::   zg1       =  2._wp        ! 
     
    114111      REAL(wp) ::   zraext_s  =  1.e+8_wp     ! extinction coefficient of radiation in the snow 
    115112      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity 
    116       REAL(wp) ::   zht_smin  =  1.e-4_wp     ! minimum snow depth 
    117113      REAL(wp) ::   ztmelt_i    ! ice melting temperature 
    118114      REAL(wp) ::   zerritmax   ! current maximal error on temperature  
     
    312308         IF( thcon_i_swi == 0 ) THEN      ! Untersteiner (1964) formula 
    313309            DO ji = kideb , kiut 
    314                ztcond_i(ji,0)        = rcdic + zbeta*s_i_b(ji,1) / MIN(-zeps,t_i_b(ji,1)-rtt) 
     310               ztcond_i(ji,0)        = rcdic + zbeta*s_i_b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt) 
    315311               ztcond_i(ji,0)        = MAX(ztcond_i(ji,0),zkimin) 
    316312            END DO 
     
    318314               DO ji = kideb , kiut 
    319315                  ztcond_i(ji,layer) = rcdic + zbeta*( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) /  & 
    320                      MIN(-2.0_wp * zeps, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) 
     316                     MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) 
    321317                  ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin) 
    322318               END DO 
     
    326322         IF( thcon_i_swi == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 
    327323            DO ji = kideb , kiut 
    328                ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( -zeps, t_i_b(ji,1)-rtt )   & 
     324               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( -epsi10, t_i_b(ji,1)-rtt )   & 
    329325                  &                   - 0.011_wp * ( t_i_b(ji,1) - rtt )   
    330326               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 
     
    333329               DO ji = kideb , kiut 
    334330                  ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) )   & 
    335                      &                                  / MIN(-2.0_wp * zeps, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)   & 
     331                     &                                  / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)   & 
    336332                     &                       - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt )   
    337333                  ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin ) 
     
    339335            END DO 
    340336            DO ji = kideb , kiut 
    341                ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-zeps,t_bo_b(ji)-rtt)   & 
     337               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt)   & 
    342338                  &                        - 0.011_wp * ( t_bo_b(ji) - rtt )   
    343339               ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 
     
    352348 
    353349            !-- Snow kappa factors 
    354             zkappa_s(ji,0)         = rcdsn / MAX(zeps,zh_s(ji)) 
    355             zkappa_s(ji,nlay_s)    = rcdsn / MAX(zeps,zh_s(ji)) 
     350            zkappa_s(ji,0)         = rcdsn / MAX(epsi10,zh_s(ji)) 
     351            zkappa_s(ji,nlay_s)    = rcdsn / MAX(epsi10,zh_s(ji)) 
    356352         END DO 
    357353 
     
    359355            DO ji = kideb , kiut 
    360356               zkappa_s(ji,layer)  = 2.0 * rcdsn / & 
    361                   MAX(zeps,2.0*zh_s(ji)) 
     357                  MAX(epsi10,2.0*zh_s(ji)) 
    362358            END DO 
    363359         END DO 
     
    367363               !-- Ice kappa factors 
    368364               zkappa_i(ji,layer)  = 2.0*ztcond_i(ji,layer)/ & 
    369                   MAX(zeps,2.0*zh_i(ji))  
    370             END DO 
    371          END DO 
    372  
    373          DO ji = kideb , kiut 
    374             zkappa_i(ji,0)        = ztcond_i(ji,0)/MAX(zeps,zh_i(ji)) 
    375             zkappa_i(ji,nlay_i)   = ztcond_i(ji,nlay_i) / MAX(zeps,zh_i(ji)) 
     365                  MAX(epsi10,2.0*zh_i(ji))  
     366            END DO 
     367         END DO 
     368 
     369         DO ji = kideb , kiut 
     370            zkappa_i(ji,0)        = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji)) 
     371            zkappa_i(ji,nlay_i)   = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji)) 
    376372            !-- Interface 
    377             zkappa_s(ji,nlay_s)   = 2.0*rcdsn*ztcond_i(ji,0)/MAX(zeps, & 
     373            zkappa_s(ji,nlay_s)   = 2.0*rcdsn*ztcond_i(ji,0)/MAX(epsi10, & 
    378374               (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji))) 
    379375            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) & 
     
    389385               ztitemp(ji,layer)   = t_i_b(ji,layer) 
    390386               zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ & 
    391                   MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),zeps) 
     387                  MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),epsi10) 
    392388               zeta_i(ji,layer)    = rdt_ice / MAX(rhoic*zspeche_i(ji,layer)*zh_i(ji), & 
    393                   zeps) 
     389                  epsi10) 
    394390            END DO 
    395391         END DO 
     
    398394            DO ji = kideb , kiut 
    399395               ztstemp(ji,layer) = t_s_b(ji,layer) 
    400                zeta_s(ji,layer)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),zeps) 
     396               zeta_s(ji,layer)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10) 
    401397            END DO 
    402398         END DO 
     
    707703      END DO  ! End of the do while iterative procedure 
    708704 
    709       IF( ln_nicep ) THEN 
     705      IF( ln_nicep .AND. lwp ) THEN 
    710706         WRITE(numout,*) ' zerritmax : ', zerritmax 
    711707         WRITE(numout,*) ' nconv     : ', nconv 
     
    732728      ! Heat conservation       ! 
    733729      !-------------------------! 
    734       IF( con_i ) THEN 
     730      IF( con_i .AND. jiindex_1d > 0 ) THEN 
    735731         DO ji = kideb, kiut 
    736732            ! Upper snow value 
Note: See TracChangeset for help on using the changeset viewer.