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

Ignore:
Timestamp:
2010-10-11T18:24:37+02:00 (14 years ago)
Author:
acc
Message:

#733 DEV_r2191_3partymerge2010. Merged in changes from DEV_r1924_nocs_latphys

File:
1 edited

Legend:

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

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