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.
stpctl.F90 in NEMO/trunk/src/OCE – NEMO

source: NEMO/trunk/src/OCE/stpctl.F90 @ 10364

Last change on this file since 10364 was 10364, checked in by acc, 5 years ago

Introduce Adaptive-Implicit vertical advection option to the trunk. This is code merged from branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP (see ticket #2042). The structure for the option is complete but is currently only successful with the flux-limited advection scheme (ln_traadv_mus). The use of this scheme with flux corrected advection schemes is not recommended until improvements to the nonoscillatory algorithm have been completed (work in progress elsewhere). The scheme is activated via a new namelist switch (ln_zad_Aimp) and is off by default.
  • Property svn:keywords set to Id
File size: 11.2 KB
RevLine 
[3]1MODULE stpctl
[1489]2   !!======================================================================
[3]3   !!                       ***  MODULE  stpctl  ***
4   !! Ocean run control :  gross check of the ocean time stepping
[1489]5   !!======================================================================
6   !! History :  OPA  ! 1991-03  (G. Madec) Original code
7   !!            6.0  ! 1992-06  (M. Imbard)
8   !!            8.0  ! 1997-06  (A.M. Treguier)
9   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
10   !!            2.0  ! 2009-07  (G. Madec)  Add statistic for time-spliting
[9019]11   !!            3.7  ! 2016-09  (G. Madec)  Remove solver
12   !!            4.0  ! 2017-04  (G. Madec)  regroup global communications
[1489]13   !!----------------------------------------------------------------------
[3]14
15   !!----------------------------------------------------------------------
16   !!   stp_ctl      : Control the run
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers variables
19   USE dom_oce         ! ocean space and time domain variables
[6140]20   USE c1d             ! 1D vertical configuration
[9210]21   USE diawri          ! Standard run outputs       (dia_wri_state routine)
[6140]22   !
[3]23   USE in_out_manager  ! I/O manager
24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
25   USE lib_mpp         ! distributed memory computing
[10364]26   USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables
[9210]27   USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy
[3]28
[9441]29   USE netcdf          ! NetCDF library
[3]30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC stp_ctl           ! routine called by step.F90
[9441]34
[10364]35   INTEGER  ::   idrun, idtime, idssh, idu, ids1, ids2, idt1, idt2, idc1, idw1, istatus
[3]36   !!----------------------------------------------------------------------
[9598]37   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1489]38   !! $Id$
[10068]39   !! Software governed by the CeCILL license (see ./LICENSE)
[1489]40   !!----------------------------------------------------------------------
[3]41CONTAINS
42
43   SUBROUTINE stp_ctl( kt, kindic )
44      !!----------------------------------------------------------------------
45      !!                    ***  ROUTINE stp_ctl  ***
46      !!                     
47      !! ** Purpose :   Control the run
48      !!
49      !! ** Method  : - Save the time step in numstp
50      !!              - Print it each 50 time steps
[9019]51      !!              - Stop the run IF problem encountered by setting indic=-3
52      !!                Problems checked: |ssh| maximum larger than 10 m
53      !!                                  |U|   maximum larger than 10 m/s
54      !!                                  negative sea surface salinity
[3]55      !!
[9019]56      !! ** Actions :   "time.step" file = last ocean time-step
57      !!                "run.stat"  file = run statistics
[9210]58      !!                nstop indicator sheared among all local domain (lk_mpp=T)
[3]59      !!----------------------------------------------------------------------
[6140]60      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
61      INTEGER, INTENT(inout) ::   kindic   ! error indicator
[1489]62      !!
[6140]63      INTEGER  ::   ji, jj, jk             ! dummy loop indices
[9019]64      INTEGER  ::   iih, ijh               ! local integers
65      INTEGER  ::   iiu, iju, iku          !   -       -
[9808]66      INTEGER  ::   iis1, ijs1, iks1       !   -       -
67      INTEGER  ::   iis2, ijs2, iks2       !   -       -
[9019]68      REAL(wp) ::   zzz                    ! local real
[9808]69      INTEGER , DIMENSION(3) ::   ilocu, ilocs1, ilocs2
[9019]70      INTEGER , DIMENSION(2) ::   iloch
[10364]71      REAL(wp), DIMENSION(9) ::   zmax
[9565]72      CHARACTER(len=20) :: clname
[3]73      !!----------------------------------------------------------------------
[6140]74      !
[3]75      IF( kt == nit000 .AND. lwp ) THEN
76         WRITE(numout,*)
77         WRITE(numout,*) 'stp_ctl : time-stepping control'
[79]78         WRITE(numout,*) '~~~~~~~'
[9019]79         !                                ! open time.step file
[1581]80         CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
[9019]81         !                                ! open run.stat file
82         CALL ctl_opn( numrun, 'run.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
[9441]83
84         IF( lwm ) THEN
[9565]85            clname = 'run.stat.nc'
86            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
87            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, idrun )
[9441]88            istatus = NF90_DEF_DIM( idrun, 'time', NF90_UNLIMITED, idtime )
89            istatus = NF90_DEF_VAR( idrun, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), idssh )
[9808]90            istatus = NF90_DEF_VAR( idrun,   'abs_u_max', NF90_DOUBLE, (/ idtime /), idu   )
91            istatus = NF90_DEF_VAR( idrun,       's_min', NF90_DOUBLE, (/ idtime /), ids1  )
92            istatus = NF90_DEF_VAR( idrun,       's_max', NF90_DOUBLE, (/ idtime /), ids2  )
[10364]93            istatus = NF90_DEF_VAR( idrun,       't_min', NF90_DOUBLE, (/ idtime /), idt1  )
94            istatus = NF90_DEF_VAR( idrun,       't_max', NF90_DOUBLE, (/ idtime /), idt2  )
95            IF( ln_zad_Aimp ) THEN
96               istatus = NF90_DEF_VAR( idrun,   'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1  )
97               istatus = NF90_DEF_VAR( idrun,       'Cu_max', NF90_DOUBLE, (/ idtime /), idc1  )
98            ENDIF
[9441]99            istatus = NF90_ENDDEF(idrun)
100         ENDIF
[10364]101         zmax(8:9) = 0._wp    ! initialise to zero in case ln_zad_Aimp option is not in use
[9441]102         
[3]103      ENDIF
[6140]104      !
[9019]105      IF(lwp) THEN                        !==  current time step  ==!   ("time.step" file)
106         WRITE ( numstp, '(1x, i8)' )   kt
107         REWIND( numstp )
[3]108      ENDIF
[6140]109      !
[9019]110      !                                   !==  test of extrema  ==!
[9023]111      IF( ll_wd ) THEN
[9210]112         zmax(1) = MAXVAL(  ABS( sshn(:,:) + ssh_ref*tmask(:,:,1) )  )        ! ssh max
[9023]113      ELSE
[9210]114         zmax(1) = MAXVAL(  ABS( sshn(:,:) )  )                               ! ssh max
[9023]115      ENDIF
[9019]116      zmax(2) = MAXVAL(  ABS( un(:,:,:) )  )                                  ! velocity max (zonal only)
117      zmax(3) = MAXVAL( -tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )   ! minus salinity max
[9808]118      zmax(4) = MAXVAL(  tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )   !       salinity max
[10364]119      zmax(5) = MAXVAL( -tsn(:,:,:,jp_tem) , mask = tmask(:,:,:) == 1._wp )   ! minus temperature max
120      zmax(6) = MAXVAL(  tsn(:,:,:,jp_tem) , mask = tmask(:,:,:) == 1._wp )   !       temperature max
121      zmax(7) = REAL( nstop , wp )                                            ! stop indicator
122      IF( ln_zad_Aimp ) THEN
123         zmax(8) = MAXVAL(  ABS( wi(:,:,:) ) , mask = wmask(:,:,:) == 1._wp ) ! implicit vertical vel. max
124         zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) !       cell Courant no. max
125      ENDIF
[1489]126      !
[9210]127      IF( lk_mpp ) THEN
[10364]128         CALL mpp_max_multiple( zmax(:), 9 )    ! max over the global domain
[9210]129         !
[10364]130         nstop = NINT( zmax(7) )                 ! nstop indicator sheared among all local domains
[9210]131      ENDIF
[1489]132      !
[9019]133      IF( MOD( kt, nwrite ) == 1 .AND. lwp ) THEN
134         WRITE(numout,*) ' ==>> time-step= ', kt, ' |ssh| max: ',   zmax(1), ' |U| max: ', zmax(2),   &
[9808]135            &                                     ' S min: '    , - zmax(3), ' S max: ', zmax(4)
[3]136      ENDIF
[6140]137      !
[9808]138      IF (  zmax(1) >   15._wp .OR.   &                    ! too large sea surface height ( > 15 m )
139         &  zmax(2) >   10._wp .OR.   &                    ! too large velocity ( > 10 m/s)
140         &  zmax(3) >=   0._wp .OR.   &                    ! negative or zero sea surface salinity
141         &  zmax(4) >= 100._wp .OR.   &                    ! too large sea surface salinity ( > 100 )
142         &  zmax(4) <    0._wp .OR.   &                    ! too large sea surface salinity (keep this line for sea-ice)
[9019]143         &  ISNAN( zmax(1) + zmax(2) + zmax(3) )  ) THEN   ! NaN encounter in the tests
144         IF( lk_mpp ) THEN
[9808]145            CALL mpp_maxloc( ABS(sshn)        , ssmask(:,:)  , zzz, iih , ijh        )
146            CALL mpp_maxloc( ABS(un)          , umask (:,:,:), zzz, iiu , iju , iku  )
147            CALL mpp_minloc( tsn(:,:,:,jp_sal), tmask (:,:,:), zzz, iis1, ijs1, iks1 )
148            CALL mpp_maxloc( tsn(:,:,:,jp_sal), tmask (:,:,:), zzz, iis2, ijs2, iks2 )
[7852]149         ELSE
[9808]150            iloch  = MINLOC( ABS( sshn(:,:)   )                               )
151            ilocu  = MAXLOC( ABS( un  (:,:,:) )                               )
152            ilocs1 = MINLOC( tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )
153            ilocs2 = MAXLOC( tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )
154            iih  = iloch (1) + nimpp - 1   ;   ijh  = iloch (2) + njmpp - 1
155            iiu  = ilocu (1) + nimpp - 1   ;   iju  = ilocu (2) + njmpp - 1   ;   iku  = ilocu (3)
156            iis1 = ilocs1(1) + nimpp - 1   ;   ijs1 = ilocs1(2) + njmpp - 1   ;   iks1 = ilocs1(3)
157            iis2 = ilocs2(1) + nimpp - 1   ;   ijs2 = ilocs2(2) + njmpp - 1   ;   iks2 = ilocs2(3)
[7852]158         ENDIF
159         IF(lwp) THEN
160            WRITE(numout,cform_err)
[9808]161            WRITE(numout,*) ' stp_ctl: |ssh| > 10 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests'
[9210]162            WRITE(numout,*) ' ======= '
[9808]163            WRITE(numout,9100) kt,   zmax(1), iih , ijh
164            WRITE(numout,9200) kt,   zmax(2), iiu , iju , iku
165            WRITE(numout,9300) kt, - zmax(3), iis1, ijs1, iks1
166            WRITE(numout,9400) kt,   zmax(4), iis2, ijs2, iks2
[7852]167            WRITE(numout,*)
[9019]168            WRITE(numout,*) '          output of last computed fields in output.abort.nc file'
[7852]169         ENDIF
170         kindic = -3
[9210]171         !
172         nstop = nstop + 1                            ! increase nstop by 1 (on all local domains)
173         CALL dia_wri_state( 'output.abort', kt )     ! create an output.abort file
174         !
[7852]175      ENDIF
[9019]1769100  FORMAT (' kt=',i8,'   |ssh| max: ',1pg11.4,', at  i j  : ',2i5)
1779200  FORMAT (' kt=',i8,'   |U|   max: ',1pg11.4,', at  i j k: ',3i5)
[9808]1789300  FORMAT (' kt=',i8,'   S     min: ',1pg11.4,', at  i j k: ',3i5)
1799400  FORMAT (' kt=',i8,'   S     max: ',1pg11.4,', at  i j k: ',3i5)
[1489]180      !
[9019]181      !                                            !==  run statistics  ==!   ("run.stat" file)
[9808]182      IF(lwp) WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3), zmax(4)
[9441]183      IF( lwm ) THEN
184         istatus = NF90_PUT_VAR( idrun, idssh, (/ zmax(1)/), (/kt/), (/1/) )
185         istatus = NF90_PUT_VAR( idrun,   idu, (/ zmax(2)/), (/kt/), (/1/) )
[9808]186         istatus = NF90_PUT_VAR( idrun,  ids1, (/-zmax(3)/), (/kt/), (/1/) )
187         istatus = NF90_PUT_VAR( idrun,  ids2, (/ zmax(4)/), (/kt/), (/1/) )
[10364]188         istatus = NF90_PUT_VAR( idrun,  idt1, (/-zmax(5)/), (/kt/), (/1/) )
189         istatus = NF90_PUT_VAR( idrun,  idt2, (/ zmax(6)/), (/kt/), (/1/) )
190         IF( ln_zad_Aimp ) THEN
191            istatus = NF90_PUT_VAR( idrun,  idw1, (/ zmax(8)/), (/kt/), (/1/) )
192            istatus = NF90_PUT_VAR( idrun,  idc1, (/ zmax(9)/), (/kt/), (/1/) )
193         ENDIF
[9441]194         IF( MOD( kt , 100 ) == 0 ) istatus = NF90_SYNC(idrun)
195         IF( kt == nitend         ) istatus = NF90_CLOSE(idrun)
196      END IF
[7852]197      !
[9808]1989500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16)
[7852]199      !
[3]200   END SUBROUTINE stp_ctl
201
202   !!======================================================================
203END MODULE stpctl
Note: See TracBrowser for help on using the repository browser.