Changeset 14357


Ignore:
Timestamp:
2021-01-28T11:45:54+01:00 (5 months ago)
Author:
cbricaud
Message:

add in dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING branch ISF top limit condition for triads (ticket #2601)

Location:
NEMO/branches/2021/dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING/src/OCE
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2021/dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING/src/OCE/LDF/ldfslp.F90

    r14329 r14357  
    148148         END_2D 
    149149      ENDIF 
    150       IF( ln_zps .AND. ln_isfcav ) THEN           ! partial steps correction at the bottom ocean level 
     150      IF( ln_zps .AND. ln_isfcav ) THEN           ! partial steps correction at the top ocean level 
    151151         DO_2D( 1, 0, 1, 0 ) 
    152152            IF( miku(ji,jj) > 1 )   zgru(ji,jj,miku(ji,jj)) = grui(ji,jj)  
     
    356356            END_2D 
    357357         ENDIF 
     358         IF( ln_zps .AND. ln_isfcav .AND. l_grad_zps) THEN           ! partial steps correction at the top ocean level 
     359            DO_2D( 1, 0, 1, 0 ) 
     360               iku  = miku(ji,jj)          ;   ikv  = mikv(ji,jj)        ! first wet ocean level (u- & v-points) 
     361               zdit = gtui(ji,jj,jp_tem)   ;   zdjt = gtvi(ji,jj,jp_tem) ! i- & j-gradient of Temperature 
     362               zdis = gtui(ji,jj,jp_sal)   ;   zdjs = gtvi(ji,jj,jp_sal) ! i- & j-gradient of Salinity 
     363               zdxrho_raw = ( - rab_b(ji+ip,jj   ,iku,jp_tem) * zdit + rab_b(ji+ip,jj   ,iku,jp_sal) * zdis ) * r1_e1u(ji,jj) 
     364               zdyrho_raw = ( - rab_b(ji   ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji   ,jj+jp,ikv,jp_sal) * zdjs ) * r1_e2v(ji,jj) 
     365               zdxrho(ji+ip,jj   ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw )   ! keep the sign 
     366               zdyrho(ji   ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 
     367            END_2D 
     368         ENDIF 
    358369         ! 
    359370      END DO 
     
    368379               zdks = 0._wp 
    369380            ENDIF 
     381            IF ( ln_isfcav ) THEN 
     382               IF( mikt(ji,jj) > 1 .AND.  jk+kp <= mikt(ji,jj) )THEN 
     383                  zdkt = 0._wp                                          ! 1st wet level gradient set to zero 
     384                  zdks = 0._wp 
     385               ENDIF 
     386            ENDIF 
     387 
    370388            zdzrho_raw = ( - rab_b(ji,jj,jk   ,jp_tem) * zdkt &  
    371389                       &   + rab_b(ji,jj,jk   ,jp_sal) * zdks & 
     
    382400      !                                       !==  intialisations to zero  ==! 
    383401      ! 
    384       wslp2  (:,:,:)     = 0._wp              ! wslp2 will be cumulated 3D field set to zero 
    385       triadi_g(:,:,1,:,:) = 0._wp   ;   triadi_g(:,:,jpk,:,:) = 0._wp   ! set surface and bottom slope to zero 
    386       triadj_g(:,:,1,:,:) = 0._wp   ;   triadj_g(:,:,jpk,:,:) = 0._wp 
    387       !!gm _iso set to zero missing 
    388       triadi  (:,:,1,:,:) = 0._wp   ;   triadj  (:,:,jpk,:,:) = 0._wp   ! set surface and bottom slope to zero 
    389       triadj  (:,:,1,:,:) = 0._wp   ;   triadj  (:,:,jpk,:,:) = 0._wp 
     402      wslp2   (:,:,:)     = 0._wp              ! wslp2 will be cumulated 3D field 
     403      triadi_g(:,:,:,:,:) = 0._wp ; triadj_g(:,:,:,:,:) = 0._wp 
     404      triadi  (:,:,:,:,:) = 0._wp ; triadj  (:,:,:,:,:) = 0._wp 
    390405 
    391406      !-------------------------------------! 
     
    445460                  ! 
    446461                  ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface 
    447                   zti_coord = znot_thru_surface * ( gdept(ji+1,jj  ,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) 
    448                   ztj_coord = znot_thru_surface * ( gdept(ji  ,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj)     ! unmasked 
     462                  zti_coord =  wmask(ji,jj,jk+kp) * znot_thru_surface * ( gdept(ji+1,jj  ,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) 
     463                  ztj_coord =  wmask(ji,jj,jk+kp) * znot_thru_surface * ( gdept(ji  ,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj)     ! unmasked 
     464 
    449465                  zti_g_raw = zti_raw - zti_coord      ! ref to geopot surfaces 
    450466                  ztj_g_raw = ztj_raw - ztj_coord 
     
    464480                  zti_g_lim =          ( zfacti   * zti_g_lim                       & 
    465481                     &      + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp)   & 
    466                      &                           * gdepw(ji+ip,jj,jk+kp,Kmm) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp) 
     482                     &                           * (gdepw(ji+ip,jj,jk+kp,Kmm)-risfdep(ji+ip,jj)) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp) 
    467483                  ztj_g_lim =          ( zfactj   * ztj_g_lim                       & 
    468484                     &      + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp)   & 
    469                      &                           * gdepw(ji,jj+jp,jk+kp,Kmm) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp) 
     485                     &                           * (gdepw(ji,jj+jp,jk+kp,Kmm)-risfdep(ji,jj+jp)) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp) 
    470486                  ! 
    471487                  triadi_g(ji+ip,jj   ,jk,1-ip,kp) = zti_g_lim 
     
    478494                  ! To do this, retain limited sx**2  in vertical flux, but divide by real slope for 13/31 terms 
    479495                  ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2 
     496                  ! 
     497                  zti_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS(zti_lim ) ), zti_lim ) 
     498                  ztj_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2, ABS(ztj_lim ) ), ztj_lim ) 
    480499                  ! 
    481500                  zti_lim  = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp)    ! remove coordinate slope => relative to coordinate surfaces 
  • NEMO/branches/2021/dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING/src/OCE/LDF/ldftra.F90

    r14329 r14357  
    248248      ENDIF 
    249249      ! 
    250       IF( ln_isfcav .AND. ln_traldf_triad )   CALL ctl_stop( ' ice shelf cavity and traldf_triad not tested' ) 
    251            ! 
    252250      IF(  nldf_tra == np_lap_i .OR. nldf_tra == np_lap_it .OR. & 
    253251         & nldf_tra == np_blp_i .OR. nldf_tra == np_blp_it  )   l_ldfslp = .TRUE.    ! slope of neutral surfaces required 
  • NEMO/branches/2021/dev_r14312_MOI-01_MOI-GRIFFIESTRIADSMOOTHING/src/OCE/TRA/traldf_triad.F90

    r14215 r14357  
    277277            !                    !==  Vertical tracer gradient at level jk and jk+1 
    278278            DO_2D( 1, 1, 1, 1 ) 
    279                zdkt3d(ji,jj,1) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1) 
     279               zdkt3d(ji,jj,1) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1) 
    280280            END_2D 
    281281            ! 
     
    284284            ELSE 
    285285               DO_2D( 1, 1, 1, 1 ) 
    286                   zdkt3d(ji,jj,0) = ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
     286                  zdkt3d(ji,jj,0) = ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) * wmask(ji,jj,jk) 
    287287               END_2D 
    288288            ENDIF 
     
    307307                        IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew 
    308308                        zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 
    309                         ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - ( zah_slp + zaei_slp) * zdxt                 * zbu * ze3wr 
     309                        ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) -                (zah_slp + zaei_slp) * zdxt   * zbu * ze3wr 
    310310                     END_2D 
    311311                  END DO 
     
    327327                        IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew 
    328328                        zftv(ji,jj   ,jk   ) = zftv(ji,jj   ,jk   ) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 
    329                         ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - ( zah_slp + zaei_slp ) * zdyt                * zbv * ze3wr 
     329                        ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) -                (zah_slp + zaei_slp) * zdyt   * zbv * ze3wr 
    330330                     END_2D 
    331331                  END DO 
     
    350350                        IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew        ! aeit(ji+ip,jj,jk)*zslope_skew 
    351351                        zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 
    352                         ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 
     352                        ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) -                (zah_slp + zaei_slp) * zdxt  * zbu * ze3wr 
    353353                     END_2D 
    354354                  END DO 
     
    369369                        zah_slp = zah * zslope_iso 
    370370                        IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew        ! aeit(ji,jj+jp,jk)*zslope_skew 
    371                         zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 
    372                         ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 
     371                        zftv(ji,jj   ,jk   ) = zftv(ji,jj   ,jk   ) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 
     372                        ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) -                (zah_slp + zaei_slp) * zdyt  * zbv * ze3wr 
    373373                     END_2D 
    374374                  END DO 
     
    388388         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz 
    389389            DO_3D( 0, 0, 1, 0, 2, jpkm1 ) 
    390                ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)   & 
     390               ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)   & 
    391391                  &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             & 
    392392                  &                            * (  pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 
     
    396396            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2 
    397397               DO_3D( 0, 0, 1, 0, 2, jpkm1 ) 
    398                   ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)             & 
     398                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)             & 
    399399                     &                            * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 
    400400               END_3D 
    401401            CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp. 
    402402               DO_3D( 0, 0, 1, 0, 2, jpkm1 ) 
    403                   ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)                      & 
     403                  ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)                      & 
    404404                     &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   & 
    405405                     &                               + akz     (ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) )   ) 
Note: See TracChangeset for help on using the changeset viewer.