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

Changeset 2163


Ignore:
Timestamp:
2010-10-06T11:30:20+02:00 (14 years ago)
Author:
sga
Message:

NEMO branch DEV_r1924_nocs_latphys. George Nurser

Included code in ldfslp.F90 to limit slopes to slpmax and
linearly decrease slopes towards surface throughout the mixed layer

Corrected minor bug for wslp2 in ldfslp.F90
Introduced rn_slpmax into namelist

Calculate streamfunctions psix_eiv,psiy_eiv

Location:
branches/DEV_r1924_nocs_latphys/NEMO/OPA_SRC
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • branches/DEV_r1924_nocs_latphys/NEMO/OPA_SRC/LDF/ldfslp.F90

    r2142 r2163  
    4141   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   alpha, beta          !: alpha,beta at T points 
    4242   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   tfw,sfw,ftu,fsu,ftv,fsv,ftud,fsud,ftvd,fsvd 
     43   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   psix_eiv 
     44   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   psiy_eiv 
    4345    
    4446   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   omlmask           ! mask of the surface mixed layer at T-pt 
    … …  
    382384          zdenr, zrhotmp, zdndt, zdddt,   &  !     "        " 
    383385          zdnds, zddds, znum, zden,       &  !     "        " 
    384           zsxe, za_sxe, zslopec, zdsloper,&  !     "        " 
     386          zslope, za_sxe, zslopec, zdsloper,&  !     "        " 
    385387          zfact, zepsln, zatempw,zatempu,zatempv, &   !     "        " 
    386           ze1ur,ze2vr,ze3wr,zdxt,zdxs,zdyt,zdys,zdzt,zdzs,zvolf 
     388          ze1ur,ze2vr,ze3wr,zdxt,zdxs,zdyt,zdys,zdzt,zdzs,zvolf,& 
     389          zr_slpmax,zdxrho,zdyrho,zabs_dzrho 
     390     REAL(wp), DIMENSION(jpi,jpj,jpk,0:1,0:1) ::   & 
     391          zsx,zsy 
     392     REAL(wp), DIMENSION(jpi,jpj,0:1,0:1) ::   & 
     393          zsx_ml_base,zsy_ml_base 
    387394     REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
    388395          zdkt,zdks 
     396     REAL(wp), DIMENSION(jpi,jpj) ::   & 
     397          zr_ml_basew 
    389398     !!---------------------------------------------------------------------- 
    390399 
    391400     ! 0. Local constant initialization 
    392401     ! -------------------------------- 
    393  
    394      zslopec=0.004 
    395      zdsloper=1000.0 
     402     zr_slpmax = 1.0_wp/slpmax 
     403 
     404     ! zslopec=0.004 
     405     ! zdsloper=1000.0 
    396406     zepsln=1e-25 
    397407 
    … …  
    495505 
    496506 
     507     ! --------------------------------------- 
     508     ! 1. Horizontal tracer gradients at T-level jk 
     509     ! --------------------------------------- 
    497510     DO jk = 1, jpkm1 
    498511        DO jj = 1, jpjm1 
    … …  
    531544        ENDIF 
    532545 
     546        ! --------------------------------------- 
     547        ! 1. Vertical tracer gradient at w-level jk 
     548        ! --------------------------------------- 
     549        DO jk = 2, jpk 
     550           zdkt(:,:,jk) = ( tb(:,:,jk-1) - tb(:,:,jk) ) * tmask(:,:,jk) 
     551           zdks(:,:,jk) = ( sb(:,:,jk-1) - sb(:,:,jk) ) * tmask(:,:,jk) 
     552        END DO 
     553 
     554        zdkt(:,:,1) = 0.0 
     555        zdks(:,:,1) = 0.0 
     556 
     557        ! --------------------------------------- 
     558        ! Depth of the w-point below ML base 
     559        ! --------------------------------------- 
     560        DO jj = 1, jpjm1 
     561           DO ji = 1, fs_jpim1 
     562              jk = nmln(ji,jj) 
     563              zr_ml_basew(ji,jj)=1.0/fsdepw(ji,jj,jk+1) 
     564           END DO 
     565        END DO 
     566 
     567 
    533568        wslp2(:,:,:)=0.0 
    534569        tfw(:,:,:) = 0.0 
    … …  
    542577        ftvd(:,:,:) = 0.0 
    543578        fsvd(:,:,:) = 0.0 
    544  
    545         DO jk = 2, jpk 
    546            ! 
    547            ! 1. Vertical tracer gradient at level jk 
    548            ! --------------------------------------- 
    549            zdkt(:,:,jk) = ( tb(:,:,jk-1) - tb(:,:,jk) ) * tmask(:,:,jk) 
    550            zdks(:,:,jk) = ( sb(:,:,jk-1) - sb(:,:,jk) ) * tmask(:,:,jk) 
    551         END DO 
    552  
    553         zdkt(:,:,1) = 0.0 
    554         zdks(:,:,1) = 0.0 
    555  
    556         ! interior (1=<jk=<jpk-1) 
     579        psix_eiv(:,:,:) = 0.0 
     580        psiy_eiv(:,:,:) = 0.0 
     581 
     582        ! ---------------------------------------------------------------------- 
     583        ! x-z plane 
     584        ! ---------------------------------------------------------------------- 
     585 
     586        ! calculate limited triad x-slopes zsx in interior (1=<jk=<jpk-1) 
    557587        DO ip=0,1 
    558588           DO kp=0,1 
    559589 
    560590              DO jk = 1, jpkm1 
    561  
    562591                 DO jj = 1, jpjm1 
    563  
    564592                    DO ji = 1, fs_jpim1   ! vector opt. 
    565593 
    … …  
    571599                       zdzt=zdkt(ji+ip,jj,jk+kp)*ze3wr 
    572600                       zdzs=zdks(ji+ip,jj,jk+kp)*ze3wr 
    573                        ! Calculate the density gradient terms 
    574                        zsxe=-(alpha(ji+ip,jj,jk)*zdxt+beta(ji+ip,jj,jk)*zdxs) / & 
    575                             (-ABS(alpha(ji+ip,jj,jk)*zdzt+beta(ji+ip,jj,jk)*zdzs)-zepsln) 
    576                        ! Check value and taper if necessary 
    577                        za_sxe = MIN(1.0,(abs(zsxe)-zslopec)*zdsloper) 
    578                        za_sxe = MAX(-1.0,za_sxe) 
    579                        zfact = 0.5*(1.0-za_sxe)*umask(ji,jj,jk)*tmask(ji+ip,jj,jk+kp) 
    580                        zsxe = zsxe*zfact 
     601                       ! Calculate the horizontal and vertical density gradient 
     602                       zdxrho = alpha(ji+ip,jj,jk)*zdxt+beta(ji+ip,jj,jk)*zdxs 
     603                       zabs_dzrho = ABS(alpha(ji+ip,jj,jk)*zdzt+beta(ji+ip,jj,jk)*zdzs)+zepsln 
     604                       ! Limit by slpmax, and mask by psi-point 
     605                       zsx(ji+ip,jj,jk,1-ip,kp) = umask(ji,jj,jk+kp) & 
     606                            & *zdxrho/MAX(zabs_dzrho,zr_slpmax*ABS(zdxrho)) 
     607                    END DO 
     608                 END DO 
     609              END DO 
     610 
     611           END DO 
     612        END DO 
     613 
     614        ! calculate slope of triad immediately below mixed-layer base 
     615        DO ip=0,1 
     616           DO kp=0,1 
     617              DO jj = 1, jpjm1 
     618                 DO ji = 1, fs_jpim1 
     619                    jk = nmln(ji,jj) 
     620                    zsx_ml_base(ji,jj,ip,kp)=zsx(ji,jj,jk+1-kp,ip,kp) 
     621                 END DO 
     622              END DO 
     623           END DO 
     624        END DO 
     625 
     626        ! Below ML use limited zsx as is 
     627        ! Inside ML replace by linearly reducing zsx_ml_base towards surface 
     628        DO ip=0,1 
     629           DO kp=0,1 
     630 
     631              DO jk = 1, jpkm1 
     632 
     633                 DO jj = 1, jpjm1 
     634 
     635                    DO ji = 1, fs_jpim1   ! vector opt. 
     636                       zfact = 1 - 1/(1 + (jk+1-kp)/nmln(ji,jj)) 
     637                       zsx(ji,jj,jk,ip,kp) = zfact*zsx(ji,jj,jk,ip,kp) + & 
     638                            & (1.0_wp-zfact)*(fsdepw(ji,jj,jk+1-kp)*zr_ml_basew(ji,jj))*zsx_ml_base(ji,jj,ip,kp)  
     639                    END DO 
     640 
     641                 END DO 
     642 
     643              END DO 
     644           END DO 
     645        END DO 
     646 
     647        ! Use zsx to calculate fluxes and save averaged slopes psix_eiv at psi-points 
     648        DO ip=0,1 
     649           DO kp=0,1 
     650 
     651              DO jk = 1, jpkm1 
     652 
     653                 DO jj = 1, jpjm1 
     654 
     655                    DO ji = 1, fs_jpim1   ! vector opt. 
     656 
     657                       ze1ur=1.0/e1u(ji,jj) 
     658                       zdxt=zdit(ji,jj,jk)*ze1ur 
     659                       zdxs=zdis(ji,jj,jk)*ze1ur 
     660 
     661                       ze3wr=1.0/fse3w(ji+ip,jj,jk+kp) 
     662                       zdzt=zdkt(ji+ip,jj,jk+kp)*ze3wr 
     663                       zdzs=zdks(ji+ip,jj,jk+kp)*ze3wr 
     664                       zslope=zsx(ji+ip,jj,jk,1-ip,kp) 
    581665 
    582666                       zvolf = 0.25_wp*e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) 
    583667 
    584                        ftu(ji,jj,jk)= ftu(ji,jj,jk)+zsxe*zdzt*zvolf*ze1ur 
    585                        fsu(ji,jj,jk)= fsu(ji,jj,jk)+zsxe*zdzs*zvolf*ze1ur 
     668                       ftu(ji,jj,jk)= ftu(ji,jj,jk)+zslope*zdzt*zvolf*ze1ur 
     669                       fsu(ji,jj,jk)= fsu(ji,jj,jk)+zslope*zdzs*zvolf*ze1ur 
    586670                       ftud(ji,jj,jk)=ftud(ji,jj,jk)+fsahtu(ji,jj,jk)*zdxt*zvolf*ze1ur 
    587671                       fsud(ji,jj,jk)=fsud(ji,jj,jk)+fsahtu(ji,jj,jk)*zdxs*zvolf*ze1ur 
    588                        tfw(ji+ip,jj,jk+kp)=tfw(ji+ip,jj,jk+kp)+(zvolf*ze3wr)*zsxe*zdxt 
    589                        sfw(ji+ip,jj,jk+kp)=sfw(ji+ip,jj,jk+kp)+(zvolf*ze3wr)*zsxe*zdxs 
     672                       tfw(ji+ip,jj,jk+kp)=tfw(ji+ip,jj,jk+kp)+(zvolf*ze3wr)*zslope*zdxt 
     673                       sfw(ji+ip,jj,jk+kp)=sfw(ji+ip,jj,jk+kp)+(zvolf*ze3wr)*zslope*zdxs 
    590674                       wslp2(ji+ip,jj,jk+kp)=wslp2(ji+ip,jj,jk+kp)+ & 
    591                             & ((zvolf*ze3wr)/(e1t(ji+ip,jj)*e2t(ji+ip,jj)))*(zsxe)**2 
     675                            & ((zvolf*ze3wr)/(e1t(ji+ip,jj)*e2t(ji+ip,jj)))*(zslope)**2 
     676                       psix_eiv(ji+ip,jj,jk+kp) = psix_eiv(ji+ip,jj,jk+kp) + 0.25_wp*zslope 
    592677 
    593678                    END DO 
    … …  
    598683        END DO 
    599684 
     685        ! ---------------------------------------------------------------------- 
     686        ! y-z plane 
     687        ! ---------------------------------------------------------------------- 
     688 
     689        ! calculate limited triad y-slopes zsy in interior (1=<jk=<jpk-1) 
    600690        DO jp=0,1 
    601691           DO kp=0,1 
    … …  
    603693              DO jk = 1, jpkm1 
    604694                 DO jj = 1, jpjm1 
    605                     DO ji = fs_2, fs_jpim1   ! vector opt. 
     695                    DO ji = 1, fs_jpim1   ! vector opt. 
    606696 
    607697                       ze2vr=1.0/e2v(ji,jj) 
    … …  
    612702                       zdzt=zdkt(ji,jj+jp,jk+kp)*ze3wr 
    613703                       zdzs=zdks(ji,jj+jp,jk+kp)*ze3wr 
    614                        ! Calculate the density gradient terms 
    615                        zsxe=-(alpha(ji,jj+jp,jk)*zdyt+beta(ji,jj+jp,jk)*zdys) / & 
    616                             (-ABS(alpha(ji,jj+jp,jk)*zdzt+beta(ji,jj+jp,jk)*zdzs)-zepsln) 
    617                        ! Check value and taper if necessary 
    618                        za_sxe = MIN(1.0,(abs(zsxe)-zslopec)*zdsloper) 
    619                        za_sxe = MAX(-1.0,za_sxe) 
    620                        zfact = 0.5*(1.0-za_sxe)*vmask(ji,jj,jk)*tmask(ji,jj+jp,jk+kp) 
    621                        zsxe = zsxe*zfact 
     704                       ! Calculate the horizontal and vertical density gradient 
     705                       zdyrho = alpha(ji,jj+jp,jk)*zdyt+beta(ji,jj+jp,jk)*zdys 
     706                       zabs_dzrho = ABS(alpha(ji,jj+jp,jk)*zdzt+beta(ji,jj+jp,jk)*zdzs)+zepsln 
     707                       ! Limit by slpmax, and mask by psi-point 
     708                       zsy(ji,jj+jp,jk,1-jp,kp) = vmask(ji,jj,jk+kp) & 
     709                            & *zdyrho/MAX(zabs_dzrho,zr_slpmax*ABS(zdyrho)) 
     710                    END DO 
     711                 END DO 
     712              END DO 
     713 
     714           END DO 
     715        END DO 
     716 
     717        ! calculate slope of triad immediately below mixed-layer base 
     718        DO jp=0,1 
     719           DO kp=0,1 
     720              DO jj = 1, jpjm1 
     721                 DO ji = 1, fs_jpim1 
     722                    jk = nmln(ji,jj) 
     723                    zsy_ml_base(ji,jj,jp,kp)=zsy(ji,jj,jk+1-kp,jp,kp) 
     724                 END DO 
     725              END DO 
     726           END DO 
     727        END DO 
     728 
     729        ! Below ML use limited zsy as is 
     730        ! Inside ML replace by linearly reducing zsy_ml_base towards surface 
     731        DO jp=0,1 
     732           DO kp=0,1 
     733 
     734              DO jk = 1, jpkm1 
     735 
     736                 DO jj = 1, jpjm1 
     737 
     738                    DO ji = 1, fs_jpim1   ! vector opt. 
     739                       zfact = 1 - 1/(1 + (jk+1-kp)/nmln(ji,jj)) 
     740                       zsy(ji,jj,jk,jp,kp) = zfact*zsy(ji,jj,jk,jp,kp) +& 
     741                            & (1.0_wp-zfact)*(fsdepw(ji,jj,jk+1-kp)*zr_ml_basew(ji,jj))*zsy_ml_base(ji,jj,jp,kp)  
     742                    END DO 
     743 
     744                 END DO 
     745 
     746              END DO 
     747           END DO 
     748        END DO 
     749 
     750        ! Use zsy to calculate fluxes and save averaged slopes psiy_eiv at psi-points 
     751        DO jp=0,1 
     752           DO kp=0,1 
     753 
     754              DO jk = 1, jpkm1 
     755 
     756                 DO jj = 1, jpjm1 
     757 
     758                    DO ji = 1, fs_jpim1   ! vector opt. 
     759 
     760                       ze2vr=1.0/e2v(ji,jj) 
     761                       zdyt=zdjt(ji,jj,jk)*ze2vr 
     762                       zdys=zdjs(ji,jj,jk)*ze2vr 
     763 
     764                       ze3wr=1.0/fse3w(ji,jj+jp,jk+kp) 
     765                       zdzt=zdkt(ji,jj+jp,jk+kp)*ze3wr 
     766                       zdzs=zdks(ji,jj+jp,jk+kp)*ze3wr 
     767                       zslope=zsy(ji,jj+jp,jk,1-jp,kp) 
    622768 
    623769                       zvolf = 0.25_wp*e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) 
    624770 
    625                        ftv(ji,jj,jk)= ftv(ji,jj,jk)+zsxe*zdzt*zvolf*ze2vr 
    626                        fsv(ji,jj,jk)= fsv(ji,jj,jk)+zsxe*zdzs*zvolf*ze2vr 
     771                       ftv(ji,jj,jk)= ftv(ji,jj,jk)+zslope*zdzt*zvolf*ze2vr 
     772                       fsv(ji,jj,jk)= fsv(ji,jj,jk)+zslope*zdzs*zvolf*ze2vr 
    627773                       ftvd(ji,jj,jk)=ftvd(ji,jj,jk)+fsahtv(ji,jj,jk)*zdyt*zvolf*ze2vr 
    628774                       fsvd(ji,jj,jk)=fsvd(ji,jj,jk)+fsahtv(ji,jj,jk)*zdys*zvolf*ze2vr 
    629                        tfw(ji,jj+jp,jk+kp)=tfw(ji,jj+jp,jk+kp)+(zvolf*ze3wr)*zsxe*zdyt 
    630                        sfw(ji,jj+jp,jk+kp)=sfw(ji,jj+jp,jk+kp)+(zvolf*ze3wr)*zsxe*zdys 
     775                       tfw(ji,jj+jp,jk+kp)=tfw(ji,jj+jp,jk+kp)+(zvolf*ze3wr)*zslope*zdyt 
     776                       sfw(ji,jj+jp,jk+kp)=sfw(ji,jj+jp,jk+kp)+(zvolf*ze3wr)*zslope*zdys 
    631777                       wslp2(ji,jj+jp,jk+kp)=wslp2(ji,jj+jp,jk+kp)+ & 
    632                             & ((zvolf*ze3wr)/(e1t(ji,jj+jp)*e2t(ji,jj+jp)))*(zsxe)**2 
     778                            & ((zvolf*ze3wr)/(e1t(ji,jj+jp)*e2t(ji,jj+jp)))*(zslope)**2 
     779                       psiy_eiv(ji,jj+jp,jk+kp) = psiy_eiv(ji,jj+jp,jk+kp) + 0.25_wp*zslope 
    633780 
    634781                    END DO 
    … …  
    654801        CALL lbc_lnk( ftvd, 'V', -1. ) 
    655802        CALL lbc_lnk( fsvd, 'V', -1. ) 
     803        CALL lbc_lnk( psix_eiv, 'U', -1. ) 
     804        CALL lbc_lnk( psiy_eiv, 'V', -1. ) 
    656805 
    657806 

  • branches/DEV_r1924_nocs_latphys/NEMO/OPA_SRC/LDF/ldftra.F90

    r2079 r2163  
    6868         &                 ln_traldf_level, ln_traldf_hor  , ln_traldf_iso,   & 
    6969         &                 ln_traldf_vis,   ln_traldf_grif ,                  & 
    70          &                 rn_aht_0       , rn_ahtb_0      , rn_aeiv_0 
     70         &                 rn_aht_0       , rn_ahtb_0      , rn_aeiv_0,       & 
     71         &                 rn_slpmax 
    7172      !!---------------------------------------------------------------------- 
    7273 
    … …  
    9293      ENDIF 
    9394 
     95      slpmax = rn_slpmax 
    9496      !                                ! convert DOCTOR namelist names into OLD names 
    9597      aht0  = rn_aht_0 

  • branches/DEV_r1924_nocs_latphys/NEMO/OPA_SRC/LDF/ldftra_oce.F90

    r2061 r2163  
    2525   REAL(wp), PUBLIC ::   rn_ahtb_0       =    0._wp  !: lateral background eddy diffusivity (m2/s) 
    2626   REAL(wp), PUBLIC ::   rn_aeiv_0       = 2000._wp  !: eddy induced velocity coefficient (m2/s) 
     27   REAL(wp), PUBLIC ::   rn_slpmax       = 0.01_wp   !: slope limit 
    2728 
    2829   REAL(wp), PUBLIC ::   aht0, ahtb0, aeiv0         !!: OLD namelist names 
    … …  
    4344   !!   'key_traldf_eiv'                              eddy induced velocity 
    4445   !!---------------------------------------------------------------------- 
    45    LOGICAL, PUBLIC, PARAMETER ::   lk_traldf_eiv   = .TRUE.   !: eddy induced velocity flag 
     46   LOGICAL, PUBLIC, PARAMETER               ::   lk_traldf_eiv   = .TRUE.   !: eddy induced velocity flag 
     47   REAL(wp), PUBLIC                         ::   slpmax                     !: slope limit  
    4648       
    4749# if defined key_traldf_c3d 

  • branches/DEV_r1924_nocs_latphys/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90

    r2060 r2163  
    178178      END DO 
    179179      ! 
     180      psix_eiv(:,:,:) = psix_eiv(:,:,:)*fsaeiu(:,:,:) 
     181      psiy_eiv(:,:,:) = psiy_eiv(:,:,:)*fsaeiv(:,:,:) 
     182 
    180183   END SUBROUTINE tra_ldf_iso_grif 
    181184 
Note: See TracChangeset for help on using the changeset viewer.