Changeset 14053 for NEMO/trunk/src/SWE/stpctl.F90

Timestamp:

2020-12-03T14:48:38+01:00 (3 years ago)

Author:

techene

Message:

#2385 added to the trunk

File:

• NEMO/trunk/src/SWE/stpctl.F90 (modified) (12 diffs)

NEMO/trunk/src/SWE/stpctl.F90

@@ -3 +3 @@
    !!                       ***  MODULE  stpctl  ***
    !! Ocean run control :  gross check of the ocean time stepping
+   !!              *** Shallow Water Equation (SWE) case ***
+   !!               ( No test on temperature and salinity )
    !!======================================================================
-   !! History :  OPA  ! 1991-03  (G. Madec) Original code
-   !!            6.0  ! 1992-06  (M. Imbard)
-   !!            8.0  ! 1997-06  (A.M. Treguier)
-   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
-   !!            2.0  ! 2009-07  (G. Madec)  Add statistic for time-spliting
-   !!            3.7  ! 2016-09  (G. Madec)  Remove solver
-   !!            4.0  ! 2017-04  (G. Madec)  regroup global communications
+   !! History :  SWE  ! 2020-09  (A. Nasser, S. Techene ) OCE/stpctl adaptated to SWE
    !!----------------------------------------------------------------------
 
@@ -21 +17 @@
    USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables
    USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy
-   !  
+   !
    USE diawri          ! Standard run outputs       (dia_wri_state routine)
    USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp         ! distributed memory computing
-   !
    USE netcdf          ! NetCDF library
+
    IMPLICIT NONE
    PRIVATE
@@ -35 +31 @@
    INTEGER                ::   nrunid   ! netcdf file id
    INTEGER, DIMENSION(2)  ::   nvarid   ! netcdf variable id
+
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -49 +47 @@
       !!
       !! ** Method  : - Save the time step in numstp
+      !!              - Print it each 50 time steps
       !!              - Stop the run IF problem encountered by setting nstop > 0
-      !!                Problems checked: negative sea surface height 
+      !!                Problems checked: e3t0+ssh minimum smaller that 0
       !!                                  |U|   maximum larger than 10 m/s 
+      !!                                  ( not for SWE : negative sea surface salinity )
       !!
       !! ** Actions :   "time.step" file = last ocean time-step
@@ -63 +63 @@
       INTEGER                         ::   idtime, istatus
       INTEGER , DIMENSION(3)          ::   iareasum, iareamin, iareamax
-      INTEGER , DIMENSION(3,2)        ::   iloc                                  ! min/max loc indices
+      INTEGER , DIMENSION(3,4)        ::   iloc                                  ! min/max loc indices
       REAL(wp)                        ::   zzz                                   ! local real 
       REAL(wp), DIMENSION(3)          ::   zmax, zmaxlocal
@@ -70 +70 @@
       CHARACTER(len=20)               ::   clname
       !!----------------------------------------------------------------------
+      IF( nstop > 0 .AND. ngrdstop > -1 )   RETURN   !   stpctl was already called by a child grid
+      !
       IF( nstop > 0 .AND. ngrdstop > -1 )   RETURN   !   stpctl was already called by a child grid
       !
@@ -109 +111 @@
       !                                   !==            test of local extrema           ==!
       !                                   !==  done by all processes at every time step  ==!
-      !
-      llmsk(   1:Nis1,:,:) = .FALSE.                                              ! exclude halos from the checked region
-      llmsk(Nie1: jpi,:,:) = .FALSE.
-      llmsk(:,   1:Njs1,:) = .FALSE.
-      llmsk(:,Nje1: jpj,:) = .FALSE.
-      !
+      zmax(1) = MINVAL( e3t_0(:,:,1)+ssh(:,:,Kmm)  )                              ! e3t_Kmm min
+      llmsk(:,:,:) = umask(:,:,:) == 1._wp
+      zmax(2) = MAXVAL(  ABS( uu(:,:,:,Kmm) ), mask = llmsk )                     ! velocity max (zonal only)
+      zmax(3) = REAL( nstop , wp )                                            ! stop indicator
+      !                                   !==               get global extrema             ==!
+      !                                   !==  done by all processes if writting run.stat  ==!
       llmsk(Nis0:Nie0,Njs0:Nje0,1) = tmask(Nis0:Nie0,Njs0:Nje0,1) == 1._wp         ! define only the inner domain
       zmax(1) = MAXVAL(     -e3t(:,:,1,Kmm) ), mask = llmsk(:,:,1) )      ! ssh max
@@ -131 +133 @@
       IF( ll_wrtruns ) THEN
          WRITE(numrun,9500) kt, zmax(1), zmax(2)
-         istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ -zmax(1)/), (/kt/), (/1/) )
-         istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/  zmax(2)/), (/kt/), (/1/) )
+         istatus = NF90_PUT_VAR( nrunid, nvarid(1), (/ zmax(1)/), (/kt/), (/1/) )
+         istatus = NF90_PUT_VAR( nrunid, nvarid(2), (/ zmax(2)/), (/kt/), (/1/) )
          IF( kt == nitend )   istatus = NF90_CLOSE(nrunid)
-      END IF
+      ENDIF
       !                                   !==               error handling               ==!
       !                                   !==  done by all processes at every time step  ==!
       !
-      IF(   zmax(1) >  0._wp           .OR.   &               ! negative sea surface height 
-         &  zmax(2) > 10._wp           .OR.   &               ! too large velocity ( > 10 m/s)
+!!SWE specific : start
+      IF(   zmax(1) <=   0._wp .OR.           &               ! negative e3t_Kmm
+         &  zmax(2) >   10._wp .OR.           &               ! too large velocity ( > 10 m/s)
          &  ISNAN( zmax(1) + zmax(2) ) .OR.   &               ! NaN encounter in the tests
          &  ABS(   zmax(1) + zmax(2) ) > HUGE(1._wp) ) THEN   ! Infinity encounter in the tests
@@ -148 +151 @@
             IF( lwm .AND. kt /= nitend )   istatus = NF90_CLOSE(nrunid)
             ! get global loc on the min/max
-            llmsk(Nis0:Nie0,Njs0:Nje0,1) = tmask(Nis0:Nie0,Njs0:Nje0,1) == 1._wp         ! define only the inner domain
-            CALL mpp_maxloc( 'stpctl',   -e3t(:,:,1,Kmm) , llmsk(:,:,1), zzz, iloc(1:2,1) )   ! mpp_maxloc ok if mask = F 
-            llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
-            CALL mpp_maxloc( 'stpctl', ABS(uu(:,:,:,Kmm)), llmsk(:,:,:), zzz, iloc(1:3,2) )
+            CALL mpp_minloc( 'stpctl', e3t_0(:,:,1) + ssh(:,:,Kmm), ssmask(:,:  ), zzz, iloc(1:2,1) )   ! mpp_maxloc ok if mask = F 
+            CALL mpp_maxloc( 'stpctl', ABS( uu(:,:,:,Kmm))        ,  umask(:,:,:), zzz, iloc(1:3,2) )
             ! find which subdomain has the max.
             iareamin(:) = jpnij+1   ;   iareamax(:) = 0   ;   iareasum(:) = 0
@@ -164 +165 @@
          ELSE                    ! find local min and max locations:
             ! if we are here, this means that the subdomain contains some oce points -> no need to test the mask used in maxloc
-            llmsk(Nis0:Nie0,Njs0:Nje0,1) = tmask(Nis0:Nie0,Njs0:Nje0,1) == 1._wp         ! define only the inner domain
-            iloc(1:2,1) = MAXLOC(   -e3t(:,:,1,Kmm) , mask = llmsk(:,:,1) )
-            llmsk(Nis0:Nie0,Njs0:Nje0,:) = umask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp        ! define only the inner domain
-            iloc(1:3,2) = MAXLOC( ABS(uu(:,:,:,Kmm)), mask = llmsk(:,:,:) )
-            DO ji = 1, 2   ! local domain indices ==> global domain indices, excluding halos
-               iloc(1:2,ji) = (/ mig0(iloc(1,ji)), mjg0(iloc(2,ji)) /)
-            END DO
+            iloc(1:2,1) = MINLOC( e3t_0(:,:,1) + ssh(:,:,Kmm), mask = ssmask(:,:  ) == 1._wp ) + (/ nimpp - 1, njmpp - 1    /)
+            iloc(1:3,2) = MAXLOC( ABS(  uu(:,:,:,       Kmm)), mask =  umask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)
             iareamin(:) = narea   ;   iareamax(:) = narea   ;   iareasum(:) = 1         ! this is local information
          ENDIF
          !
-         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests'
-         CALL wrt_line( ctmp2, kt, '|e3t| min', -zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) )
-         CALL wrt_line( ctmp3, kt, '|U|   max',  zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) )
+         WRITE(ctmp1,*) ' stp_ctl:  e3t0+ssh < 0 m  or  |U| > 10 m/s  or  NaN encounter in the tests'
+         CALL wrt_line( ctmp2, kt, 'e3t0+ssh min',  zmax(1), iloc(:,1), iareasum(1), iareamin(1), iareamax(1) )
+         CALL wrt_line( ctmp3, kt, '|U|   max'   ,  zmax(2), iloc(:,2), iareasum(2), iareamin(2), iareamax(2) )
          IF( Agrif_Root() ) THEN
             WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort* files'
@@ -194 +190 @@
          !
       ENDIF
+!!SWE specific : end
       !
       IF( nstop > 0 ) THEN                                                  ! an error was detected and we did not abort yet...
@@ -200 +197 @@
       ENDIF
       !
-9500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16)
+9500  FORMAT(' it :', i8, '      e3t_min: ', D23.16, ' |U|_max: ', D23.16)
       !
    END SUBROUTINE stp_ctl