Changeset 3666 for branches/2012/dev_LOCEAN_UKMO_CMCC_INGV_2012/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90

Timestamp:

2012-11-26T15:22:04+01:00 (11 years ago)

Author:

cetlod

Message:

commit the changes resulting for the merged branches, see ticket #1025

File:

• branches/2012/dev_LOCEAN_UKMO_CMCC_INGV_2012/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90 (modified) (12 diffs)

branches/2012/dev_LOCEAN_UKMO_CMCC_INGV_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) include the upstream where needed
    !!----------------------------------------------------------------------
 
@@ -18 +19 @@
    USE trdmod_oce      ! tracers trends 
    USE trdtra      ! tracers trends 
+   USE eosbn2          ! equation of state 
    USE in_out_manager  ! I/O manager
    USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
    USE trabbl          ! tracers: bottom boundary layer
+   USE sbcrnf          ! river runoffs
    USE lib_mpp         ! distribued memory computing
    USE lbclnk          ! ocean lateral boundary condition (or mpp link) 
@@ -27 +30 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE eosbn2          ! equation of state
+   USE sbcrnf          ! river runoffs
 
    IMPLICIT NONE
@@ -35 +40 @@
    LOGICAL  :: l_trd       ! flag to compute trends
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits
+   !                                                             !  and in closed seas (orca 2 and 4 configurations)
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -78 +85 @@
       REAL(wp) ::   ztra, zbtr, zdt, zalpha   !   -      -
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zslpx, zslpy
+      INTEGER  ::   ierr
+      REAL(wp) ::   zice                             ! temporary scalars
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztfreez
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zind
       !!----------------------------------------------------------------------
       !
@@ -83 +94 @@
       !
       CALL wrk_alloc( jpi, jpj, jpk, zslpx, zslpy )
+      CALL wrk_alloc( jpi, jpj, ztfreez )
+      CALL wrk_alloc( jpi, jpj, jpk, zind )
       !
 
@@ -89 +102 @@
          IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~'
+         IF(lwp) WRITE(numout,*)
+         !
+         !
+         IF (.not. ALLOCATED(upsmsk))THEN
+             ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
+             IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_muscl: unable to allocate array')
+         ENDIF
+         !
+         upsmsk(:,:) = 0._wp                             ! not upstream by default
          !
          l_trd = .FALSE.
@@ -94 +116 @@
       ENDIF
 
+      !
+      ! Upstream / centered scheme indicator
+      ! ------------------------------------
+      ztfreez(:,:) = tfreez( tsn(:,:,1,jp_sal) )
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !                                        ! below ice covered area (if tn < "freezing"+0.1 )
+               IF( tsn(ji,jj,jk,jp_tem) <= ztfreez(ji,jj) + 0.1_wp ) THEN   ;   zice = 1.e0
+               ELSE                                                         ;   zice = 0.e0
+               ENDIF
+               zind(ji,jj,jk) = MAX (   &
+                  rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows)
+                  upsmsk(ji,jj)               ,      &  ! some of some straits
+                  zice                               &  ! below ice covered area (if tn < "freezing"+0.1 )
+                  &                  ) * tmask(ji,jj,jk)
+               zind(ji,jj,jk) = 1 - zind(ji,jj,jk)
+            END DO
+         END DO
+      END DO
       !                                                     ! ===========
       DO jn = 1, kjpt                                       ! tracer loop
@@ -148 +190 @@
                   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 )
                   !
@@ -155 +197 @@
                   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
@@ -230 +272 @@
                   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 
@@ -255 +297 @@
       !
       CALL wrk_dealloc( jpi, jpj, jpk, zslpx, zslpy )
+      CALL wrk_dealloc( jpi, jpj, ztfreez )
+      CALL wrk_dealloc( jpi, jpj, jpk, zind )
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_adv_muscl')