Changeset 3680 for branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90

Timestamp:

2012-11-27T15:42:24+01:00 (11 years ago)

Author:

rblod

Message:

First commit of the final branch for 2012 (future nemo_3_5), see ticket #1028

File:

• branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90 (modified) (8 diffs)

branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90

@@ -8 +8 @@
    !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
    !!            3.2  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
+   !!            3.4  !  2012-06  (P. Oddo, M. Vichi) include the upstream where needed
    !!----------------------------------------------------------------------
 
@@ -28 +29 @@
    USE timing         ! Timing
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE eosbn2          ! equation of state
+   USE sbcrnf          ! river runoffs
 
    IMPLICIT NONE
@@ -34 +37 @@
    PUBLIC   tra_adv_muscl  ! routine called by step.F90
 
-   LOGICAL  :: l_trd       ! flag to compute trends
-
+   LOGICAL  :: l_trd                        ! flag to compute trends
+   LOGICAL, PUBLIC  :: ln_traadv_msc_ups= .FALSE.   ! use upstream scheme within muscl
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits
+   !                                                             !  and in closed seas (orca 2 and 4 configurations)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zind         !: mixed upstream/centered index
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -79 +86 @@
       REAL(wp) ::   ztra, zbtr, zdt, zalpha   !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zslpx, zslpy
+      INTEGER  ::   ierr
       !!----------------------------------------------------------------------
       !
@@ -89 +97 @@
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
+         IF(lwp) WRITE(numout,*) '        : xed up-stream            ' , ln_traadv_msc_ups
          IF(lwp) WRITE(numout,*) '~~~~~~~'
+         IF(lwp) WRITE(numout,*)
+         !
+         !
+         IF(ln_traadv_msc_ups) THEN
+           IF (.not. ALLOCATED(upsmsk))THEN
+               ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
+               IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_muscl: unable to allocate upsmsk array')
+           ENDIF
+           upsmsk(:,:) = 0._wp                             ! not upstream by default
+         ENDIF
+
+         IF (.not. ALLOCATED(zind))THEN
+             ALLOCATE( zind(jpi,jpj,jpk), STAT=ierr )
+             IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_muscl: unable to allocate zind array')
+         ENDIF
+         !
          !
          l_trd = .FALSE.
          IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
-      ENDIF
-
+
+      !
+      ! Upstream / centered scheme indicator
+      ! ------------------------------------
+         zind(:,:,:) = 1._wp                             ! set equal to 0 where up-stream is needed
+
+         IF(ln_traadv_msc_ups) THEN
+           DO jk = 1, jpk
+              DO jj = 1, jpj
+                 DO ji = 1, jpi
+                    zind(ji,jj,jk) = 1  - MAX (           &
+                       rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows)
+                       upsmsk(ji,jj) ) * tmask(ji,jj,jk)     ! some of some straits
+                 END DO
+              END DO
+           END DO
+         ENDIF 
+      !
+      ENDIF ! end kit000
       !                                                     ! ===========
       DO jn = 1, kjpt                                       ! tracer loop
@@ -149 +191 @@
                   zalpha = 0.5 - z0u
                   zu  = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
-                  zzwx = ptb(ji+1,jj,jk,jn) + zu * zslpx(ji+1,jj,jk)
-                  zzwy = ptb(ji  ,jj,jk,jn) + zu * zslpx(ji  ,jj,jk)
+                  zzwx = ptb(ji+1,jj,jk,jn) + zind(ji,jj,jk) * (zu * zslpx(ji+1,jj,jk))
+                  zzwy = ptb(ji  ,jj,jk,jn) + zind(ji,jj,jk) * (zu * zslpx(ji  ,jj,jk))
                   zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
                   !
@@ -156 +198 @@
                   zalpha = 0.5 - z0v
                   zv  = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
-                  zzwx = ptb(ji,jj+1,jk,jn) + zv * zslpy(ji,jj+1,jk)
-                  zzwy = ptb(ji,jj  ,jk,jn) + zv * zslpy(ji,jj  ,jk) 
+                  zzwx = ptb(ji,jj+1,jk,jn) + zind(ji,jj,jk) * (zv * zslpy(ji,jj+1,jk))
+                  zzwy = ptb(ji,jj  ,jk,jn) + zind(ji,jj,jk) * (zv * zslpy(ji,jj  ,jk))
                   zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
                END DO
@@ -231 +273 @@
                   zalpha = 0.5 + z0w
                   zw  = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr 
-                  zzwx = ptb(ji,jj,jk+1,jn) + zw * zslpx(ji,jj,jk+1)
-                  zzwy = ptb(ji,jj,jk  ,jn) + zw * zslpx(ji,jj,jk  )
+                  zzwx = ptb(ji,jj,jk+1,jn) + zind(ji,jj,jk) * (zw * zslpx(ji,jj,jk+1))
+                  zzwy = ptb(ji,jj,jk  ,jn) + zind(ji,jj,jk) * (zw * zslpx(ji,jj,jk  ))
                   zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
                END DO