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 @ 10370

Last change on this file since 10370 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
Line 
1MODULE stpctl
2   !!======================================================================
3   !!                       ***  MODULE  stpctl  ***
4   !! Ocean run control :  gross check of the ocean time stepping
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
11   !!            3.7  ! 2016-09  (G. Madec)  Remove solver
12   !!            4.0  ! 2017-04  (G. Madec)  regroup global communications
13   !!----------------------------------------------------------------------
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
20   USE c1d             ! 1D vertical configuration
21   USE diawri          ! Standard run outputs       (dia_wri_state routine)
22   !
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
26   USE zdf_oce ,  ONLY : ln_zad_Aimp       ! ocean vertical physics variables
27   USE wet_dry,   ONLY : ll_wd, ssh_ref    ! reference depth for negative bathy
28
29   USE netcdf          ! NetCDF library
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC stp_ctl           ! routine called by step.F90
34
35   INTEGER  ::   idrun, idtime, idssh, idu, ids1, ids2, idt1, idt2, idc1, idw1, istatus
36   !!----------------------------------------------------------------------
37   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
38   !! $Id$
39   !! Software governed by the CeCILL license (see ./LICENSE)
40   !!----------------------------------------------------------------------
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
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
55      !!
56      !! ** Actions :   "time.step" file = last ocean time-step
57      !!                "run.stat"  file = run statistics
58      !!                nstop indicator sheared among all local domain (lk_mpp=T)
59      !!----------------------------------------------------------------------
60      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
61      INTEGER, INTENT(inout) ::   kindic   ! error indicator
62      !!
63      INTEGER  ::   ji, jj, jk             ! dummy loop indices
64      INTEGER  ::   iih, ijh               ! local integers
65      INTEGER  ::   iiu, iju, iku          !   -       -
66      INTEGER  ::   iis1, ijs1, iks1       !   -       -
67      INTEGER  ::   iis2, ijs2, iks2       !   -       -
68      REAL(wp) ::   zzz                    ! local real
69      INTEGER , DIMENSION(3) ::   ilocu, ilocs1, ilocs2
70      INTEGER , DIMENSION(2) ::   iloch
71      REAL(wp), DIMENSION(9) ::   zmax
72      CHARACTER(len=20) :: clname
73      !!----------------------------------------------------------------------
74      !
75      IF( kt == nit000 .AND. lwp ) THEN
76         WRITE(numout,*)
77         WRITE(numout,*) 'stp_ctl : time-stepping control'
78         WRITE(numout,*) '~~~~~~~'
79         !                                ! open time.step file
80         CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
81         !                                ! open run.stat file
82         CALL ctl_opn( numrun, 'run.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
83
84         IF( lwm ) THEN
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 )
88            istatus = NF90_DEF_DIM( idrun, 'time', NF90_UNLIMITED, idtime )
89            istatus = NF90_DEF_VAR( idrun, 'abs_ssh_max', NF90_DOUBLE, (/ idtime /), idssh )
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  )
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
99            istatus = NF90_ENDDEF(idrun)
100         ENDIF
101         zmax(8:9) = 0._wp    ! initialise to zero in case ln_zad_Aimp option is not in use
102         
103      ENDIF
104      !
105      IF(lwp) THEN                        !==  current time step  ==!   ("time.step" file)
106         WRITE ( numstp, '(1x, i8)' )   kt
107         REWIND( numstp )
108      ENDIF
109      !
110      !                                   !==  test of extrema  ==!
111      IF( ll_wd ) THEN
112         zmax(1) = MAXVAL(  ABS( sshn(:,:) + ssh_ref*tmask(:,:,1) )  )        ! ssh max
113      ELSE
114         zmax(1) = MAXVAL(  ABS( sshn(:,:) )  )                               ! ssh max
115      ENDIF
116      zmax(2) = MAXVAL(  ABS( un(:,:,:) )  )                                  ! velocity max (zonal only)
117      zmax(3) = MAXVAL( -tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )   ! minus salinity max
118      zmax(4) = MAXVAL(  tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )   !       salinity max
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
126      !
127      IF( lk_mpp ) THEN
128         CALL mpp_max_multiple( zmax(:), 9 )    ! max over the global domain
129         !
130         nstop = NINT( zmax(7) )                 ! nstop indicator sheared among all local domains
131      ENDIF
132      !
133      IF( MOD( kt, nwrite ) == 1 .AND. lwp ) THEN
134         WRITE(numout,*) ' ==>> time-step= ', kt, ' |ssh| max: ',   zmax(1), ' |U| max: ', zmax(2),   &
135            &                                     ' S min: '    , - zmax(3), ' S max: ', zmax(4)
136      ENDIF
137      !
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)
143         &  ISNAN( zmax(1) + zmax(2) + zmax(3) )  ) THEN   ! NaN encounter in the tests
144         IF( lk_mpp ) THEN
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 )
149         ELSE
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)
158         ENDIF
159         IF(lwp) THEN
160            WRITE(numout,cform_err)
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'
162            WRITE(numout,*) ' ======= '
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
167            WRITE(numout,*)
168            WRITE(numout,*) '          output of last computed fields in output.abort.nc file'
169         ENDIF
170         kindic = -3
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         !
175      ENDIF
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)
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)
180      !
181      !                                            !==  run statistics  ==!   ("run.stat" file)
182      IF(lwp) WRITE(numrun,9500) kt, zmax(1), zmax(2), -zmax(3), zmax(4)
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/) )
186         istatus = NF90_PUT_VAR( idrun,  ids1, (/-zmax(3)/), (/kt/), (/1/) )
187         istatus = NF90_PUT_VAR( idrun,  ids2, (/ zmax(4)/), (/kt/), (/1/) )
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
194         IF( MOD( kt , 100 ) == 0 ) istatus = NF90_SYNC(idrun)
195         IF( kt == nitend         ) istatus = NF90_CLOSE(idrun)
196      END IF
197      !
1989500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16)
199      !
200   END SUBROUTINE stp_ctl
201
202   !!======================================================================
203END MODULE stpctl
Note: See TracBrowser for help on using the repository browser.