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 14283 – NEMO

Changeset 14283


Ignore:
Timestamp:
2021-01-08T18:33:46+01:00 (3 years ago)
Author:
agn
Message:

correct various bugs in shear code

zd_cubic def in working code
correct bugs found by checking reproducibility
correct bugs found in email of 26 Nov 20
correct loop index limit in zdf_osm_osbl_state
Ensure salinity in ztsm_midu does not go to zero in zdf_osm_zmld_horizontal_gradients
Corrected bug Simon found in PE calcuklation

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90

    r14281 r14283  
    715715      DO jj = 2, jpjm1 
    716716        DO ji = 2, jpim1 
    717           IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) + jp_ext(ji,jj) >= mbkt(ji,jj) .or. ibld(ji,jj)-imld(jj,ji) == 1 ) lpyc(ji,jj) = .FALSE.  
     717          IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) + jp_ext(ji,jj) >= mbkt(ji,jj) .or. ibld(ji,jj)-imld(ji,jj) == 1 ) lpyc(ji,jj) = .FALSE.  
    718718        END DO 
    719719      END DO         
     
    915915             IF ( j_ddh(ji,jj) == 0 ) THEN 
    916916! Place holding code. Parametrization needs checking for these conditions. 
    917                zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird 
     917               zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) ))**pthird 
    918918               zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) 
    919919               zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) 
    920920             ELSE 
    921                zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird 
     921               zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) ))**pthird 
    922922               zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj) 
    923923               zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj) 
     
    928928             DO jk = imld(ji,jj), ibld(ji,jj) 
    929929                zznd_pyc = - ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj) 
    930                 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 
     930                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045 * ( ztau_sc_u(ji,jj)**2 ) * zuw_bse * & 
     931                     & ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 
    931932             END DO 
    932933             zvw_max = 0.7 * ff_t(ji,jj) * ( zustke(ji,jj) * dstokes(ji,jj) + 0.75 * zustar(ji,jj) * zhml(ji,jj) ) 
     
    935936             DO jk = imld(ji,jj), ibld(ji,jj) 
    936937               zznd_pyc = -( gdepw_n(ji,jj,jk) -zhbl(ji,jj) ) / zdh(ji,jj) 
    937                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 
     938               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045 * ( ztau_sc_u(ji,jj)**2 ) * zvw_bse *  & 
     939                    & ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk) 
    938940             END DO 
    939941           ENDIF  ! lpyc 
     
    10671069            IF ( .not. lconv(ji,jj) ) THEN 
    10681070               DO jk = 2, ibld(ji,jj) 
    1069                   znd = ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about 
     1071                  znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about 
    10701072                  IF ( znd >= 0.0 ) THEN 
    10711073                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) 
     
    10861088             IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN 
    10871089                DO jk= 2, ibld(ji,jj) 
    1088                    znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj) 
    10891090                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 
    10901091                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) 
     
    14071408                 zbeta_v_sc(ji,jj) = 1.0 -  2.0 * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln ) 
    14081409               ELSE 
    1409                  zbeta_d_sc = 1.0 
    1410                  zbeta_v_sc = 1.0 
     1410                 zbeta_d_sc(ji,jj) = 1.0 
     1411                 zbeta_v_sc(ji,jj) = 1.0 
    14111412               ENDIF 
    14121413             ELSE 
     
    14481449                   zb_cubic = -1.75 * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj) 
    14491450                   zd_cubic = ( 0.5 * zvisml_sc(ji,jj) * zdh(ji,jj) / zhml(ji,jj) - 0.85 * zvispyc_s_sc(ji,jj)  )  / MAX(zvispyc_n_sc(ji,jj), 1.e-8) 
    1450                    zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zd_cubic ) 
     1451                   zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zb_cubic ) 
    14511452                   zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic 
    14521453                   DO jk = imld(ji,jj) , ibld(ji,jj) 
     
    16041605 
    16051606     DO jj = 2, jpjm1 
    1606        DO ji = 2, jpjm1 
     1607       DO ji = 2, jpim1 
    16071608         IF ( lshear(ji,jj) ) THEN 
    16081609           IF ( lconv(ji,jj) ) THEN 
     
    17791780                zpe_mle_layer = 0._wp 
    17801781                zpe_mle_ref = 0._wp 
     1782                zthermal = rab_n(ji,jj,1,jp_tem) 
     1783                zbeta    = rab_n(ji,jj,1,jp_sal) 
     1784 
    17811785                DO jk = ibld(ji,jj), mld_prof(ji,jj) 
    17821786                  zbuoy = grav * ( zthermal * tsn(ji,jj,jk,jp_tem) - zbeta * tsn(ji,jj,jk,jp_sal) ) 
     
    18161820                     lpyc(ji,jj) = .FALSE. 
    18171821                   ELSE 
    1818                       lpyc = .TRUE. 
     1822                      lpyc(ji,jj) = .TRUE. 
    18191823                   ENDIF 
    18201824                 ELSE 
     
    25032507      END DO 
    25042508 
     2509      ! ensure salinity > 0 in unset values so EOS doesn't give FP error with fpe0 on 
     2510      ztsm_midu(:,jpj,jp_sal) = 10. 
     2511      ztsm_midv(jpi,:,jp_sal) = 10. 
     2512 
    25052513      CALL eos_rab(ztsm_midu, zmld_midu, zabu) 
    25062514      CALL eos_rab(ztsm_midv, zmld_midv, zabv) 
Note: See TracChangeset for help on using the changeset viewer.