Ignore:
Timestamp:
2020-02-25T16:29:34+01:00 (12 months ago)
Author:
jcastill
Message:

First implementation of the branch - compiling after merge

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/DYN/dynnxt.F90

    r11715 r12453  
    3333   USE dynadv         ! dynamics: vector invariant versus flux form 
    3434   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme 
     35   USE dynspg 
    3536   USE domvvl         ! variable volume 
    3637   USE bdy_oce   , ONLY: ln_bdy 
     
    178179            vn(:,:,jk) = va(:,:,jk) 
    179180         END DO 
     181         ! limit velocities  
     182         IF (ln_ulimit) THEN  
     183            call dyn_limit_velocity (kt)  
     184         ENDIF 
    180185         IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n 
    181186!!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a   
     
    210215               END DO 
    211216            END DO 
     217            ! limit velocities  
     218            IF (ln_ulimit) THEN  
     219               call dyn_limit_velocity (kt)  
     220            ENDIF 
    212221            !                             ! ================! 
    213222         ELSE                             ! Variable volume ! 
     
    275284                  END DO 
    276285               END DO 
     286               ! limit velocities  
     287               IF (ln_ulimit) THEN  
     288                  call dyn_limit_velocity (kt)  
     289               ENDIF 
    277290               ! 
    278291            ELSE                          ! Asselin filter applied on thickness weighted velocity 
     
    302315                  END DO 
    303316               END DO 
     317               ! limit velocities  
     318               IF (ln_ulimit) THEN  
     319                  call dyn_limit_velocity (kt)  
     320               ENDIF 
    304321               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor 
    305322               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 
     
    377394   END SUBROUTINE dyn_nxt 
    378395 
     396   SUBROUTINE dyn_limit_velocity (kt)  
     397      ! limits maxming vlaues of un and vn by volume courant number  
     398      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index  
     399      !  
     400      INTEGER  ::   ji, jj, jk   ! dummy loop indices  
     401      REAL(wp) ::   zzu,zplim,zmlim,isp,ism,zcn,ze3e1,zzcn,zcnn,idivp,idivm  
     402  
     403      ! limit fluxes  
     404      zcn =cn_ulimit !0.9 ! maximum velocity inverse courant number  
     405      zcnn = cnn_ulimit !0.54 ! how much to reduce cn by in divergen flow  
     406  
     407      DO jk = 1, jpkm1  
     408         DO jj = 1, jpjm1  
     409            DO ji = 1, jpim1  
     410               ! U direction  
     411               zzu = un(ji,jj,jk)  
     412               ze3e1 = e3u_n(ji  ,jj,jk) * e2u(ji  ,jj)  
     413               ! ips is 1 if flow is positive othersize zero  
     414               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp )  
     415               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp )  
     416               ! idev is 1 if divergent flow otherwise zero  
     417               idivp = -isp * 0.5 * (sign(1.0_wp, un(ji-1,jj,jk)) - 1.0_wp )  
     418               idivm =  ism * 0.5 * (sign(1.0_wp, un(ji+1,jj,jk)) + 1.0_wp )  
     419               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp  
     420               zzcn = zzcn * zcn  
     421               zplim =  zzcn * (e3t_n(ji  ,jj,jk) * e1t(ji  ,jj) * e2t(ji  ,jj)) / &  
     422                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk)  
     423               zmlim = -zzcn * (e3t_n(ji+1,jj,jk) * e1t(ji+1,jj) * e2t(ji+1,jj)) / &  
     424                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk)  
     425               ! limit currents  
     426               un(ji,jj,jk) = min ( zzu,zplim) * isp + max(zzu,zmlim) *ism  
     427               ! V  direction  
     428               zzu = vn(ji,jj,jk)  
     429               ze3e1 = e3v_n(ji  ,jj,jk) * e1v(ji  ,jj)  
     430               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp )  
     431               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp )  
     432               ! idev is 1 if divergent flow otherwise zero  
     433               idivp = -isp * 0.5 * (sign(1.0_wp, vn(ji,jj-1,jk)) - 1.0_wp )  
     434               idivm =  ism * 0.5 * (sign(1.0_wp, vn(ji,jj+1,jk)) + 1.0_wp )  
     435               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp  
     436               zzcn = zzcn * zcn  
     437               zplim =  zzcn * (e3t_n(ji,jj  ,jk) * e1t(ji,jj  ) * e2t(ji,jj  )) / &  
     438                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk)  
     439               zmlim = -zzcn * (e3t_n(ji,jj+1,jk) * e1t(ji,jj+1) * e2t(ji,jj+1)) / &  
     440                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk)  
     441               vn(ji,jj,jk) = min ( zzu,zplim) * isp + max(zzu,zmlim) *ism  
     442            ENDDO  
     443         ENDDO  
     444      ENDDO  
     445  
     446    END SUBROUTINE dyn_limit_velocity 
     447 
    379448   !!========================================================================= 
    380449END MODULE dynnxt 
Note: See TracChangeset for help on using the changeset viewer.