source: trunk/NEMOGCM/NEMO/OPA_SRC/stpctl.F90 @ 7852

Last change on this file since 7852 was 7852, checked in by cetlod, 4 years ago

trunk:clear stop with output.abort created after numerical explosion of ssh, see shaconemo ticket #66

  • Property svn:keywords set to Id
File size: 7.8 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   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   stp_ctl      : Control the run
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers variables
17   USE dom_oce         ! ocean space and time domain variables
18   USE c1d             ! 1D vertical configuration
19   !
20   USE in_out_manager  ! I/O manager
21   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
22   USE lib_mpp         ! distributed memory computing
23   USE lib_fortran     ! Fortran routines library
24
25   IMPLICIT NONE
26   PRIVATE
27
28   PUBLIC stp_ctl           ! routine called by step.F90
29   !!----------------------------------------------------------------------
30   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
31   !! $Id$
32   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
33   !!----------------------------------------------------------------------
34CONTAINS
35
36   SUBROUTINE stp_ctl( kt, kindic )
37      !!----------------------------------------------------------------------
38      !!                    ***  ROUTINE stp_ctl  ***
39      !!                     
40      !! ** Purpose :   Control the run
41      !!
42      !! ** Method  : - Save the time step in numstp
43      !!              - Print it each 50 time steps
44      !!              - Stop the run IF problem ( indic < 0 )
45      !!
46      !! ** Actions :   'time.step' file containing the last ocean time-step
47      !!               
48      !!----------------------------------------------------------------------
49      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
50      INTEGER, INTENT(inout) ::   kindic   ! error indicator
51      !!
52      INTEGER  ::   ji, jj, jk             ! dummy loop indices
53      INTEGER  ::   ii, ij, ik             ! local integers
54      REAL(wp) ::   zumax, zsmin, zssh2, zsshmax    ! local scalars
55      INTEGER, DIMENSION(3) ::   ilocu     !
56      INTEGER, DIMENSION(2) ::   ilocs     !
57      !!----------------------------------------------------------------------
58      !
59      IF( kt == nit000 .AND. lwp ) THEN
60         WRITE(numout,*)
61         WRITE(numout,*) 'stp_ctl : time-stepping control'
62         WRITE(numout,*) '~~~~~~~'
63         ! open time.step file
64         CALL ctl_opn( numstp, 'time.step', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
65      ENDIF
66      !
67      IF(lwp) WRITE ( numstp, '(1x, i8)' )   kt      !* save the current time step in numstp
68      IF(lwp) REWIND( numstp )                       !  --------------------------
69      !
70      !                                              !* Test maximum of velocity (zonal only)
71      !                                              !  ------------------------
72      !! zumax = MAXVAL( ABS( un(:,:,:) ) )                ! slower than the following loop on NEC SX5
73      zumax = 0.e0
74      DO jk = 1, jpk
75         DO jj = 1, jpj
76            DO ji = 1, jpi
77               zumax = MAX(zumax,ABS(un(ji,jj,jk)))
78          END DO
79        END DO
80      END DO       
81      IF( lk_mpp )   CALL mpp_max( zumax )                 ! max over the global domain
82      !
83      IF( MOD( kt, nwrite ) == 1 .AND. lwp )   WRITE(numout,*) ' ==>> time-step= ',kt,' abs(U) max: ', zumax
84      !
85      IF( zumax > 20.e0 ) THEN
86         IF( lk_mpp ) THEN
87            CALL mpp_maxloc(ABS(un),umask,zumax,ii,ij,ik)
88         ELSE
89            ilocu = MAXLOC( ABS( un(:,:,:) ) )
90            ii = ilocu(1) + nimpp - 1
91            ij = ilocu(2) + njmpp - 1
92            ik = ilocu(3)
93         ENDIF
94         IF(lwp) THEN
95            WRITE(numout,cform_err)
96            WRITE(numout,*) ' stpctl: the zonal velocity is larger than 20 m/s'
97            WRITE(numout,*) ' ====== '
98            WRITE(numout,9400) kt, zumax, ii, ij, ik
99            WRITE(numout,*)
100            WRITE(numout,*) '          output of last fields in numwso'
101         ENDIF
102         kindic = -3
103      ENDIF
1049400  FORMAT (' kt=',i6,' max abs(U): ',1pg11.4,', i j k: ',3i5)
105      !
106      !                                              !* Test minimum of salinity
107      !                                              !  ------------------------
108      !! zsmin = MINVAL( tsn(:,:,1,jp_sal), mask = tmask(:,:,1) == 1.e0 )  slower than the following loop on NEC SX5
109      zsmin = 100._wp
110      DO jj = 2, jpjm1
111         DO ji = 1, jpi
112            IF( tmask(ji,jj,1) == 1) zsmin = MIN(zsmin,tsn(ji,jj,1,jp_sal))
113         END DO
114      END DO
115      IF( lk_mpp )   CALL mpp_min( zsmin )                ! min over the global domain
116      !
117      IF( MOD( kt, nwrite ) == 1 .AND. lwp )   WRITE(numout,*) ' ==>> time-step= ',kt,' SSS min:', zsmin
118      !
119      IF( zsmin < 0.) THEN
120         IF (lk_mpp) THEN
121            CALL mpp_minloc ( tsn(:,:,1,jp_sal),tmask(:,:,1), zsmin, ii,ij )
122         ELSE
123            ilocs = MINLOC( tsn(:,:,1,jp_sal), mask = tmask(:,:,1) == 1.e0 )
124            ii = ilocs(1) + nimpp - 1
125            ij = ilocs(2) + njmpp - 1
126         ENDIF
127         !
128         IF(lwp) THEN
129            WRITE(numout,cform_err)
130            WRITE(numout,*) 'stp_ctl : NEGATIVE sea surface salinity'
131            WRITE(numout,*) '======= '
132            WRITE(numout,9500) kt, zsmin, ii, ij
133            WRITE(numout,*)
134            WRITE(numout,*) '          output of last fields in numwso'
135         ENDIF
136         kindic = -3
137      ENDIF
1389500  FORMAT (' kt=',i6,' min SSS: ',1pg11.4,', i j: ',2i5)
139      !
140      !
141      IF( lk_c1d )  RETURN          ! No log file in case of 1D vertical configuration
142
143      ! log file (ssh statistics)
144      ! --------                                   !* ssh statistics (and others...)
145      IF( kt == nit000 .AND. lwp ) THEN   ! open ssh statistics file (put in solver.stat file)
146         CALL ctl_opn( numsol, 'solver.stat', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
147      ENDIF
148      !
149      zsshmax = 0.e0
150      DO jj = 1, jpj
151         DO ji = 1, jpi
152            IF( tmask(ji,jj,1) == 1) zsshmax = MAX( zsshmax, ABS(sshn(ji,jj)) )
153         END DO
154      END DO
155      IF( lk_mpp )   CALL mpp_max( zsshmax )                ! min over the global domain
156      !
157      IF( MOD( kt, nwrite ) == 1 .AND. lwp )   WRITE(numout,*) ' ==>> time-step= ',kt,' ssh max:', zsshmax
158      !
159      IF( zsshmax > 10.e0 ) THEN
160         IF (lk_mpp) THEN
161            CALL mpp_maxloc( ABS(sshn(:,:)),tmask(:,:,1),zsshmax,ii,ij)
162         ELSE
163            ilocs = MAXLOC( ABS(sshn(:,:)) )
164            ii = ilocs(1) + nimpp - 1
165            ij = ilocs(2) + njmpp - 1
166         ENDIF
167         !
168         IF(lwp) THEN
169            WRITE(numout,cform_err)
170            WRITE(numout,*) 'stp_ctl : the ssh is larger than 10m'
171            WRITE(numout,*) '======= '
172            WRITE(numout,9600) kt, zsshmax, ii, ij
173            WRITE(numout,*)
174            WRITE(numout,*) '          output of last fields in numwso'
175         ENDIF
176         kindic = -3
177      ENDIF
1789600  FORMAT (' kt=',i6,' max ssh: ',1pg11.4,', i j: ',2i5)
179      !
180      zssh2 = glob_sum( sshn(:,:) * sshn(:,:) )
181      !
182      IF(lwp) WRITE(numsol,9700) kt, zssh2, zumax, zsmin      ! ssh statistics
183      !
1849700  FORMAT(' it :', i8, ' ssh2: ', d23.16, ' Umax: ',d23.16,' Smin: ',d23.16)
185      !
186   END SUBROUTINE stp_ctl
187
188   !!======================================================================
189END MODULE stpctl
Note: See TracBrowser for help on using the repository browser.