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 2236 for branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/LDF/ldfslp.F90 – NEMO

Ignore:
Timestamp:
2010-10-12T20:49:32+02:00 (14 years ago)
Author:
cetlod
Message:

First guess of NEMO_v3.3

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/LDF/ldfslp.F90

    r2104 r2236  
    2424   USE phycst          ! physical constants 
    2525   USE zdfmxl          ! mixed layer depth 
     26   USE eosbn2 
    2627   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2728   USE in_out_manager  ! I/O manager 
     
    3334   PUBLIC   ldf_slp         ! routine called by step.F90 
    3435   PUBLIC   ldf_slp_init    ! routine called by opa.F90 
     36   PUBLIC ldf_slp_grif   !              " 
    3537 
    3638   LOGICAL , PUBLIC, PARAMETER              ::   lk_ldfslp = .TRUE.   !: slopes flag 
    3739   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   uslp, wslpi          !: i_slope at U- and W-points 
    3840   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   vslp, wslpj          !: j-slope at V- and W-points 
     41   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   wslp2                !: wslp**2 from Griffies quarter cells 
     42   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   alpha, beta          !: alpha,beta at T points 
     43   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   tfw,sfw,ftu,fsu,ftv,fsv,ftud,fsud,ftvd,fsvd 
     44   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   psix_eiv 
     45   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   psiy_eiv 
    3946    
    4047   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   omlmask           ! mask of the surface mixed layer at T-pt 
     
    4451   !! * Substitutions 
    4552#  include "domzgr_substitute.h90" 
     53#  include "ldftra_substitute.h90" 
     54#  include "ldfeiv_substitute.h90" 
    4655#  include "vectopt_loop_substitute.h90" 
    4756   !!---------------------------------------------------------------------- 
    4857   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)  
    4958   !! $Id$ 
    50    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     59   !! Software governed by the CeCILL licence  (NEMOGCM/License_CeCILL.txt) 
    5160   !!---------------------------------------------------------------------- 
    5261 
     
    304313      END DO 
    305314          
    306  
    307315      ! III.  Specific grid points      
    308316      ! ===========================  
     
    339347      ! 
    340348   END SUBROUTINE ldf_slp 
     349 
     350   SUBROUTINE ldf_slp_grif ( kt ) 
     351     !!---------------------------------------------------------------------- 
     352     !!                 ***  ROUTINE ldf_slp_grif  *** 
     353     !! 
     354     !! ** Purpose :   Compute the squared slopes of neutral surfaces (slope 
     355     !!      of iso-pycnal surfaces referenced locally) ('key_traldfiso') 
     356     !!      at W-points using the Griffies quarter-cells.  Also calculates 
     357     !!      alpha and beta at T-points for use in the Griffies isopycnal 
     358     !!      scheme. 
     359     !! 
     360     !! ** Method  :   The slope in the i-direction is computed at U- and 
     361     !!      W-points (uslp, wslpi) and the slope in the j-direction is 
     362     !!      computed at V- and W-points (vslp, wslpj). 
     363     !! 
     364     !! ** Action : - alpha, beta 
     365     !!               wslp2 squared slope of neutral surfaces at w-points. 
     366     !! 
     367     !! History : 
     368     !!   9.0  !  06-10  (C. Harris)  New subroutine 
     369     !!---------------------------------------------------------------------- 
     370     !! * Modules used 
     371     USE oce            , zdit  => ua,  &  ! use ua as workspace 
     372          zdis  => va,  &  ! use va as workspace 
     373          zdjt  => ta,  &  ! use ta as workspace 
     374          zdjs  => sa      ! use sa as workspace 
     375     !! * Arguments 
     376     INTEGER, INTENT( in ) ::   kt         ! ocean time-step index 
     377 
     378     !! * Local declarations 
     379     INTEGER  ::   ji, jj, jk, ip, jp, kp  ! dummy loop indices 
     380     INTEGER  ::   iku, ikv                ! temporary integer 
     381     REAL(wp) ::   & 
     382          zt, zs, zh, zt2, zsp5, zp1t1,   &  ! temporary scalars 
     383          zdenr, zrhotmp, zdndt, zdddt,   &  !     "        " 
     384          zdnds, zddds, znum, zden,       &  !     "        " 
     385          zslope, za_sxe, zslopec, zdsloper,&  !     "        " 
     386          zfact, zepsln, zatempw,zatempu,zatempv, &   !     "        " 
     387          ze1ur,ze2vr,ze3wr,zdxt,zdxs,zdyt,zdys,zdzt,zdzs,zvolf,& 
     388          zr_slpmax,zdxrho,zdyrho,zabs_dzrho 
     389     REAL(wp), DIMENSION(jpi,jpj,jpk,0:1,0:1) ::   & 
     390          zsx,zsy 
     391     REAL(wp), DIMENSION(jpi,jpj,0:1,0:1) ::   & 
     392          zsx_ml_base,zsy_ml_base 
     393     REAL(wp), DIMENSION(jpi,jpj,jpk) ::   & 
     394          zdkt,zdks 
     395     REAL(wp), DIMENSION(jpi,jpj) ::   & 
     396          zr_ml_basew 
     397     !!---------------------------------------------------------------------- 
     398 
     399     ! 0. Local constant initialization 
     400     ! -------------------------------- 
     401     zr_slpmax = 1.0_wp/slpmax 
     402 
     403     ! zslopec=0.004 
     404     ! zdsloper=1000.0 
     405     zepsln=1e-25 
     406 
     407     SELECT CASE ( nn_eos ) 
     408 
     409     CASE ( 0 )               ! Jackett and McDougall (1994) formulation 
     410 
     411        !                                                ! =============== 
     412        DO jk = 1, jpk                                   ! Horizontal slab 
     413           !                                             ! =============== 
     414           DO jj = 1, jpjm1 
     415              DO ji = 1, fs_jpim1 
     416                 zt = tb(ji,jj,jk)          ! potential temperature 
     417                 zs = sb(ji,jj,jk) - 35.0   ! salinity anomaly (s-35) 
     418                 zh = fsdept(ji,jj,jk)      ! depth in meters 
     419 
     420                 beta(ji,jj,jk) = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt      & 
     421                      &                            - 0.301985e-05 ) * zt      & 
     422                      &   + 0.785567e-03                                      & 
     423                      &   + (     0.515032e-08 * zs                           & 
     424                      &         + 0.788212e-08 * zt - 0.356603e-06 ) * zs     & 
     425                      &   +(  (   0.121551e-17 * zh                           & 
     426                      &         - 0.602281e-15 * zs                           & 
     427                      &         - 0.175379e-14 * zt + 0.176621e-12 ) * zh     & 
     428                      &                             + 0.408195e-10   * zs     & 
     429                      &     + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt     & 
     430                      &                             - 0.121555e-07 ) * zh 
     431 
     432                 alpha(ji,jj,jk) = - beta(ji,jj,jk) *                       & 
     433                      &     (((( - 0.255019e-07 * zt + 0.298357e-05 ) * zt   & 
     434                      &                               - 0.203814e-03 ) * zt   & 
     435                      &                               + 0.170907e-01 ) * zt   & 
     436                      &   + 0.665157e-01                                      & 
     437                      &   +     ( - 0.678662e-05 * zs                         & 
     438                      &           - 0.846960e-04 * zt + 0.378110e-02 ) * zs   & 
     439                      &   +   ( ( - 0.302285e-13 * zh                         & 
     440                      &           - 0.251520e-11 * zs                         & 
     441                      &           + 0.512857e-12 * zt * zt           ) * zh   & 
     442                      &           - 0.164759e-06 * zs                         & 
     443                      &        +(   0.791325e-08 * zt - 0.933746e-06 ) * zt   & 
     444                      &                               + 0.380374e-04 ) * zh) 
     445 
     446              ENDDO 
     447           ENDDO 
     448        ENDDO 
     449 
     450     CASE ( 1 ) 
     451 
     452        alpha(:,:,:)=-rn_alpha 
     453        beta(:,:,:)=0.0 
     454 
     455     CASE ( 2 ) 
     456 
     457        alpha(:,:,:)=-rn_alpha 
     458        beta (:,:,:)=rn_beta 
     459 
     460     CASE ( 3 ) 
     461 
     462        DO jk =1, jpk 
     463           DO jj = 1, jpjm1 
     464              DO ji = 1, fs_jpim1 
     465                 zt = tb(ji,jj,jk) 
     466                 zs = sb(ji,jj,jk) 
     467                 zh = fsdept(ji,jj,jk) 
     468                 zt2 = zt**2 
     469                 zsp5 = sqrt(ABS(zs)) 
     470                 zp1t1=zh*zt 
     471                 znum=9.99843699e+02+zt*(7.35212840e+00+zt*(-5.45928211e-02+3.98476704e-04*zt)) & 
     472                      +zs*(2.96938239e+00-7.23268813e-03*zt+2.12382341e-03*zs)  & 
     473                      +zh*(1.04004591e-02+1.03970529e-07*zt2+5.18761880e-06*zs+ & 
     474                      zh*(-3.24041825e-08-1.23869360e-11*zt2)) 
     475                 zden=1.00000000e+00+zt*(7.28606739e-03+zt*(-4.60835542e-05+zt*(3.68390573e-07+zt*1.80809186e-10))) & 
     476                      +zs*(2.14691708e-03+zt*(-9.27062484e-06-1.78343643e-10*zt2)+zsp5*(4.76534122e-06+1.63410736e-09*zt2)) & 
     477                      + zh*(5.30848875e-06+zh*zt*(-3.03175128e-16*zt2-1.27934137e-17*zh)) 
     478                 zdenr=1.0/zden 
     479                 zrhotmp=znum*zdenr 
     480                 zdndt=7.35212840e+00+zt*(-1.091856422e-01+1.195430112e-03*zt)-7.23268813e-03*zs & 
     481                      +zp1t1*(2.07941058e-07-2.4773872e-11*zh) 
     482                 zdddt=7.28606739e-03+zt*(-9.21671084e-05+zt*(1.105171719e-06+7.23236744e-10*zt))  & 
     483                      +zs*(-9.27062484e-06+zt*(-5.35030929e-10*zt+3.26821472e-09*zsp5))  & 
     484                      +zh*zh*(-9.09525384e-16*zt2-1.27934137e-17*zh) 
     485                 zdnds=2.96938239e+00-7.23268813e-03*zt+2*2.12382341e-03*zs+5.18761880e-06*zh 
     486                 zddds=2.14691708e-03+zt*(-9.27062484e-06-1.78343643e-10*zt2)+zsp5*(7.14801183e-06+2.45116104e-09*zt2) 
     487                 alpha(ji,jj,jk)=(zdndt-zrhotmp*zdddt)*zdenr 
     488                 beta(ji,jj,jk)=zdenr*(zdnds-zrhotmp*zddds) 
     489 
     490              END DO 
     491           END DO 
     492        END DO 
     493 
     494     CASE DEFAULT 
     495 
     496        IF(lwp) WRITE(numout,cform_err) 
     497        IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos 
     498        nstop = nstop + 1 
     499 
     500     END SELECT 
     501 
     502     CALL lbc_lnk( alpha, 'T', 1. ) 
     503     CALL lbc_lnk( beta, 'T', 1. ) 
     504 
     505 
     506     ! --------------------------------------- 
     507     ! 1. Horizontal tracer gradients at T-level jk 
     508     ! --------------------------------------- 
     509     DO jk = 1, jpkm1 
     510        DO jj = 1, jpjm1 
     511           DO ji = 1, fs_jpim1   ! vector opt. 
     512              ! i-gradient of T and S at jj 
     513              zdit (ji,jj,jk) = ( tb(ji+1,jj,jk)-tb(ji,jj,jk) ) * umask(ji,jj,jk) 
     514              zdis (ji,jj,jk) = ( sb(ji+1,jj,jk)-sb(ji,jj,jk) ) * umask(ji,jj,jk) 
     515              ! j-gradient of T and S at jj 
     516              zdjt (ji,jj,jk) = ( tb(ji,jj+1,jk)-tb(ji,jj,jk) ) * vmask(ji,jj,jk) 
     517              zdjs (ji,jj,jk) = ( sb(ji,jj+1,jk)-sb(ji,jj,jk) ) * vmask(ji,jj,jk) 
     518           END DO 
     519        END DO 
     520     END DO 
     521 
     522     IF( ln_zps ) THEN      ! partial steps correction at the last level 
     523# if defined key_vectopt_loop   &&   ! defined key_mpp_omp 
     524     jj = 1 
     525        DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
     526# else 
     527           DO jj = 1, jpjm1 
     528              DO ji = 1, jpim1 
     529# endif 
     530                 ! last ocean level 
     531                 iku  = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1 
     532                 ikv  = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1 
     533                 ! i-gradient of T and S 
     534                 zdit (ji,jj,iku) = gtu(ji,jj) 
     535                 zdis (ji,jj,iku) = gsu(ji,jj) 
     536                 ! j-gradient of T and S 
     537                 zdjt (ji,jj,ikv) = gtv(ji,jj) 
     538                 zdjs (ji,jj,ikv) = gsv(ji,jj) 
     539# if ! defined key_vectopt_loop   ||   defined key_mpp_omp 
     540              END DO 
     541# endif 
     542           END DO 
     543        ENDIF 
     544 
     545        ! --------------------------------------- 
     546        ! 1. Vertical tracer gradient at w-level jk 
     547        ! --------------------------------------- 
     548        DO jk = 2, jpk 
     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        ! --------------------------------------- 
     557        ! Depth of the w-point below ML base 
     558        ! --------------------------------------- 
     559        DO jj = 1, jpj 
     560           DO ji = 1, jpi 
     561              jk = nmln(ji,jj) 
     562              zr_ml_basew(ji,jj)=1.0/fsdepw(ji,jj,jk+1) 
     563           END DO 
     564        END DO 
     565 
     566 
     567        wslp2(:,:,:)=0.0 
     568        tfw(:,:,:) = 0.0 
     569        sfw(:,:,:) = 0.0 
     570        ftu(:,:,:) = 0.0 
     571        fsu(:,:,:) = 0.0 
     572        ftv(:,:,:) = 0.0 
     573        fsv(:,:,:) = 0.0 
     574        ftud(:,:,:) = 0.0 
     575        fsud(:,:,:) = 0.0 
     576        ftvd(:,:,:) = 0.0 
     577        fsvd(:,:,:) = 0.0 
     578        psix_eiv(:,:,:) = 0.0 
     579        psiy_eiv(:,:,:) = 0.0 
     580 
     581        ! ---------------------------------------------------------------------- 
     582        ! x-z plane 
     583        ! ---------------------------------------------------------------------- 
     584 
     585        ! calculate limited triad x-slopes zsx in interior (1=<jk=<jpk-1) 
     586        DO ip=0,1 
     587           DO kp=0,1 
     588 
     589              DO jk = 1, jpkm1 
     590                 DO jj = 1, jpjm1 
     591                    DO ji = 1, fs_jpim1   ! vector opt. 
     592 
     593                       ze1ur=1.0/e1u(ji,jj) 
     594                       zdxt=zdit(ji,jj,jk)*ze1ur 
     595                       zdxs=zdis(ji,jj,jk)*ze1ur 
     596 
     597                       ze3wr=1.0/fse3w(ji+ip,jj,jk+kp) 
     598                       zdzt=zdkt(ji+ip,jj,jk+kp)*ze3wr 
     599                       zdzs=zdks(ji+ip,jj,jk+kp)*ze3wr 
     600                       ! Calculate the horizontal and vertical density gradient 
     601                       zdxrho = alpha(ji+ip,jj,jk)*zdxt+beta(ji+ip,jj,jk)*zdxs 
     602                       zabs_dzrho = ABS(alpha(ji+ip,jj,jk)*zdzt+beta(ji+ip,jj,jk)*zdzs)+zepsln 
     603                       ! Limit by slpmax, and mask by psi-point 
     604                       zsx(ji+ip,jj,jk,1-ip,kp) = umask(ji,jj,jk+kp) & 
     605                            & *zdxrho/MAX(zabs_dzrho,zr_slpmax*ABS(zdxrho)) 
     606                    END DO 
     607                 END DO 
     608              END DO 
     609 
     610           END DO 
     611        END DO 
     612 
     613        ! calculate slope of triad immediately below mixed-layer base 
     614        DO ip=0,1 
     615           DO kp=0,1 
     616              DO jj = 1, jpjm1 
     617                 DO ji = 1, fs_jpim1 
     618                    jk = nmln(ji+ip,jj) 
     619                    zsx_ml_base(ji+ip,jj,1-ip,kp)=zsx(ji+ip,jj,jk+1-kp,1-ip,kp) 
     620                 END DO 
     621              END DO 
     622           END DO 
     623        END DO 
     624 
     625        ! Below ML use limited zsx as is 
     626        ! Inside ML replace by linearly reducing zsx_ml_base towards surface 
     627        DO ip=0,1 
     628           DO kp=0,1 
     629 
     630              DO jk = 1, jpkm1 
     631 
     632                 DO jj = 1, jpjm1 
     633 
     634                    DO ji = 1, fs_jpim1   ! vector opt. 
     635                       ! k index of uppermost point(s) of triad is jk+kp-1 
     636                       ! must be .ge. nmln(ji,jj) for zfact=1. 
     637                       !                   otherwise  zfact=0. 
     638                       zfact = 1 - 1/(1 + (jk+kp-1)/nmln(ji+ip,jj)) 
     639                       zsx(ji+ip,jj,jk,1-ip,kp) = zfact*zsx(ji+ip,jj,jk,1-ip,kp) + & 
     640                            & (1.0_wp-zfact)*(fsdepw(ji+ip,jj,jk+kp)*zr_ml_basew(ji+ip,jj))*zsx_ml_base(ji+ip,jj,1-ip,kp)  
     641                    END DO 
     642 
     643                 END DO 
     644 
     645              END DO 
     646           END DO 
     647        END DO 
     648 
     649        ! Use zsx to calculate fluxes and save averaged slopes psix_eiv at psi-points 
     650        DO ip=0,1 
     651           DO kp=0,1 
     652 
     653              DO jk = 1, jpkm1 
     654 
     655                 DO jj = 1, jpjm1 
     656 
     657                    DO ji = 1, fs_jpim1   ! vector opt. 
     658 
     659                       ze1ur=1.0/e1u(ji,jj) 
     660                       zdxt=zdit(ji,jj,jk)*ze1ur 
     661                       zdxs=zdis(ji,jj,jk)*ze1ur 
     662 
     663                       ze3wr=1.0/fse3w(ji+ip,jj,jk+kp) 
     664                       zdzt=zdkt(ji+ip,jj,jk+kp)*ze3wr 
     665                       zdzs=zdks(ji+ip,jj,jk+kp)*ze3wr 
     666                       zslope=zsx(ji+ip,jj,jk,1-ip,kp) 
     667 
     668                       zvolf = 0.25_wp*e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) 
     669 
     670                       ftu(ji,jj,jk)= ftu(ji,jj,jk)+zslope*zdzt*zvolf*ze1ur 
     671                       fsu(ji,jj,jk)= fsu(ji,jj,jk)+zslope*zdzs*zvolf*ze1ur 
     672                       ftud(ji,jj,jk)=ftud(ji,jj,jk)+fsahtu(ji,jj,jk)*zdxt*zvolf*ze1ur 
     673                       fsud(ji,jj,jk)=fsud(ji,jj,jk)+fsahtu(ji,jj,jk)*zdxs*zvolf*ze1ur 
     674                       tfw(ji+ip,jj,jk+kp)=tfw(ji+ip,jj,jk+kp)+(zvolf*ze3wr)*zslope*zdxt 
     675                       sfw(ji+ip,jj,jk+kp)=sfw(ji+ip,jj,jk+kp)+(zvolf*ze3wr)*zslope*zdxs 
     676                       wslp2(ji+ip,jj,jk+kp)=wslp2(ji+ip,jj,jk+kp)+ & 
     677                            & ((zvolf*ze3wr)/(e1t(ji+ip,jj)*e2t(ji+ip,jj)))*(zslope)**2 
     678                       psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp*zslope 
     679 
     680                    END DO 
     681                 END DO 
     682 
     683              END DO 
     684           END DO 
     685        END DO 
     686 
     687        ! ---------------------------------------------------------------------- 
     688        ! y-z plane 
     689        ! ---------------------------------------------------------------------- 
     690 
     691        ! calculate limited triad y-slopes zsy in interior (1=<jk=<jpk-1) 
     692        DO jp=0,1 
     693           DO kp=0,1 
     694 
     695              DO jk = 1, jpkm1 
     696                 DO jj = 1, jpjm1 
     697                    DO ji = 1, fs_jpim1   ! vector opt. 
     698 
     699                       ze2vr=1.0/e2v(ji,jj) 
     700                       zdyt=zdjt(ji,jj,jk)*ze2vr 
     701                       zdys=zdjs(ji,jj,jk)*ze2vr 
     702 
     703                       ze3wr=1.0/fse3w(ji,jj+jp,jk+kp) 
     704                       zdzt=zdkt(ji,jj+jp,jk+kp)*ze3wr 
     705                       zdzs=zdks(ji,jj+jp,jk+kp)*ze3wr 
     706                       ! Calculate the horizontal and vertical density gradient 
     707                       zdyrho = alpha(ji,jj+jp,jk)*zdyt+beta(ji,jj+jp,jk)*zdys 
     708                       zabs_dzrho = ABS(alpha(ji,jj+jp,jk)*zdzt+beta(ji,jj+jp,jk)*zdzs)+zepsln 
     709                       ! Limit by slpmax, and mask by psi-point 
     710                       zsy(ji,jj+jp,jk,1-jp,kp) = vmask(ji,jj,jk+kp) & 
     711                            & *zdyrho/MAX(zabs_dzrho,zr_slpmax*ABS(zdyrho)) 
     712                    END DO 
     713                 END DO 
     714              END DO 
     715 
     716           END DO 
     717        END DO 
     718 
     719        ! calculate slope of triad immediately below mixed-layer base 
     720        DO jp=0,1 
     721           DO kp=0,1 
     722              DO jj = 1, jpjm1 
     723                 DO ji = 1, fs_jpim1 
     724                    jk = nmln(ji,jj+jp) 
     725                    zsy_ml_base(ji,jj+jp,1-jp,kp)=zsy(ji,jj+jp,jk+1-kp,1-jp,kp) 
     726                 END DO 
     727              END DO 
     728           END DO 
     729        END DO 
     730 
     731        ! Below ML use limited zsy as is 
     732        ! Inside ML replace by linearly reducing zsy_ml_base towards surface 
     733        DO jp=0,1 
     734           DO kp=0,1 
     735 
     736              DO jk = 1, jpkm1 
     737 
     738                 DO jj = 1, jpjm1 
     739 
     740                    DO ji = 1, fs_jpim1   ! vector opt. 
     741                       ! k index of uppermost point(s) of triad is jk+kp-1 
     742                       ! must be .ge. nmln(ji,jj) for zfact=1. 
     743                       !                   otherwise  zfact=0. 
     744                       zfact = 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)) 
     745                       zsy(ji,jj+jp,jk,1-jp,kp) = zfact*zsy(ji,jj+jp,jk,1-jp,kp) + & 
     746                            & (1.0_wp-zfact)*(fsdepw(ji,jj+jp,jk+kp)*zr_ml_basew(ji,jj+jp))*zsy_ml_base(ji,jj+jp,1-jp,kp)  
     747                    END DO 
     748 
     749                 END DO 
     750 
     751              END DO 
     752           END DO 
     753        END DO 
     754 
     755        ! Use zsy to calculate fluxes and save averaged slopes psiy_eiv at psi-points 
     756        DO jp=0,1 
     757           DO kp=0,1 
     758 
     759              DO jk = 1, jpkm1 
     760 
     761                 DO jj = 1, jpjm1 
     762 
     763                    DO ji = 1, fs_jpim1   ! vector opt. 
     764 
     765                       ze2vr=1.0/e2v(ji,jj) 
     766                       zdyt=zdjt(ji,jj,jk)*ze2vr 
     767                       zdys=zdjs(ji,jj,jk)*ze2vr 
     768 
     769                       ze3wr=1.0/fse3w(ji,jj+jp,jk+kp) 
     770                       zdzt=zdkt(ji,jj+jp,jk+kp)*ze3wr 
     771                       zdzs=zdks(ji,jj+jp,jk+kp)*ze3wr 
     772                       zslope=zsy(ji,jj+jp,jk,1-jp,kp) 
     773 
     774                       zvolf = 0.25_wp*e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) 
     775 
     776                       ftv(ji,jj,jk)= ftv(ji,jj,jk)+zslope*zdzt*zvolf*ze2vr 
     777                       fsv(ji,jj,jk)= fsv(ji,jj,jk)+zslope*zdzs*zvolf*ze2vr 
     778                       ftvd(ji,jj,jk)=ftvd(ji,jj,jk)+fsahtv(ji,jj,jk)*zdyt*zvolf*ze2vr 
     779                       fsvd(ji,jj,jk)=fsvd(ji,jj,jk)+fsahtv(ji,jj,jk)*zdys*zvolf*ze2vr 
     780                       tfw(ji,jj+jp,jk+kp)=tfw(ji,jj+jp,jk+kp)+(zvolf*ze3wr)*zslope*zdyt 
     781                       sfw(ji,jj+jp,jk+kp)=sfw(ji,jj+jp,jk+kp)+(zvolf*ze3wr)*zslope*zdys 
     782                       wslp2(ji,jj+jp,jk+kp)=wslp2(ji,jj+jp,jk+kp)+ & 
     783                            & ((zvolf*ze3wr)/(e1t(ji,jj+jp)*e2t(ji,jj+jp)))*(zslope)**2 
     784                       psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp*zslope 
     785 
     786                    END DO 
     787                 END DO 
     788 
     789              END DO 
     790           END DO 
     791        END DO 
     792 
     793        tfw(:,:,1)=0.0 
     794        sfw(:,:,1)=0.0 
     795        wslp2(:,:,1)=0.0 
     796 
     797        CALL lbc_lnk( wslp2, 'W', 1. ) 
     798        CALL lbc_lnk( tfw, 'W', 1. ) 
     799        CALL lbc_lnk( sfw, 'W', 1. ) 
     800        CALL lbc_lnk( ftu, 'U', -1. ) 
     801        CALL lbc_lnk( fsu, 'U', -1. ) 
     802        CALL lbc_lnk( ftv, 'V', -1. ) 
     803        CALL lbc_lnk( fsv, 'V', -1. ) 
     804        CALL lbc_lnk( ftud, 'U', -1. ) 
     805        CALL lbc_lnk( fsud, 'U', -1. ) 
     806        CALL lbc_lnk( ftvd, 'V', -1. ) 
     807        CALL lbc_lnk( fsvd, 'V', -1. ) 
     808        CALL lbc_lnk( psix_eiv, 'U', -1. ) 
     809        CALL lbc_lnk( psiy_eiv, 'V', -1. ) 
     810 
     811 
     812      END SUBROUTINE ldf_slp_grif 
    341813 
    342814 
Note: See TracChangeset for help on using the changeset viewer.