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.
dia25h.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/dia25h.F90 @ 10965

Last change on this file since 10965 was 10965, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : DIA and stpctl.F90. Just testing in ORCA1 so far.
  • Property svn:keywords set to Id
File size: 12.3 KB
Line 
1MODULE dia25h 
2   !!======================================================================
3   !!                       ***  MODULE  diaharm  ***
4   !! Harmonic analysis of tidal constituents
5   !!======================================================================
6   !! History :  3.6  !  2014  (E O'Dea)  Original code
7   !!----------------------------------------------------------------------
8   USE oce             ! ocean dynamics and tracers variables
9   USE dom_oce         ! ocean space and time domain
10   USE zdf_oce         ! ocean vertical physics   
11   USE zdfgls   , ONLY : hmxl_n
12   !
13   USE in_out_manager  ! I/O units
14   USE iom             ! I/0 library
15   USE wet_dry
16
17   IMPLICIT NONE
18   PRIVATE
19
20   PUBLIC   dia_25h_init               ! routine called by nemogcm.F90
21   PUBLIC   dia_25h                    ! routine called by diawri.F90
22
23   LOGICAL, PUBLIC ::   ln_dia25h      !:  25h mean output
24
25   ! variables for calculating 25-hourly means
26   INTEGER , SAVE ::   cnt_25h           ! Counter for 25 hour means
27   REAL(wp), SAVE ::   r1_25 = 0.04_wp   ! =1/25
28   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   tn_25h  , sn_25h
29   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   sshn_25h 
30   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   un_25h  , vn_25h  , wn_25h
31   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avt_25h , avm_25h
32   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   en_25h  , rmxln_25h
33
34   !!----------------------------------------------------------------------
35   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
36   !! $Id$
37   !! Software governed by the CeCILL license (see ./LICENSE)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE dia_25h_init( Kbb )
42      !!---------------------------------------------------------------------------
43      !!                  ***  ROUTINE dia_25h_init  ***
44      !!     
45      !! ** Purpose: Initialization of 25h mean namelist
46      !!       
47      !! ** Method : Read namelist
48      !!---------------------------------------------------------------------------
49      INTEGER, INTENT(in) :: Kbb       ! Time level index
50      !
51      INTEGER ::   ios                 ! Local integer output status for namelist read
52      INTEGER ::   ierror              ! Local integer for memory allocation
53      !
54      NAMELIST/nam_dia25h/ ln_dia25h
55      !!----------------------------------------------------------------------
56      !
57      REWIND ( numnam_ref )              ! Read Namelist nam_dia25h in reference namelist : 25hour mean diagnostics
58      READ   ( numnam_ref, nam_dia25h, IOSTAT=ios, ERR= 901 )
59901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_dia25h in reference namelist', lwp )
60      REWIND( numnam_cfg )              ! Namelist nam_dia25h in configuration namelist  25hour diagnostics
61      READ  ( numnam_cfg, nam_dia25h, IOSTAT = ios, ERR = 902 )
62902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_dia25h in configuration namelist', lwp )
63      IF(lwm) WRITE ( numond, nam_dia25h )
64
65      IF(lwp) THEN                   ! Control print
66         WRITE(numout,*)
67         WRITE(numout,*) 'dia_25h_init : Output 25 hour mean diagnostics'
68         WRITE(numout,*) '~~~~~~~~~~~~'
69         WRITE(numout,*) '   Namelist nam_dia25h : set 25h outputs '
70         WRITE(numout,*) '      Switch for 25h diagnostics (T) or not (F)  ln_dia25h  = ', ln_dia25h
71      ENDIF
72      IF( .NOT. ln_dia25h )   RETURN
73      ! ------------------- !
74      ! 1 - Allocate memory !
75      ! ------------------- !
76      !                                ! ocean arrays
77      ALLOCATE( tn_25h (jpi,jpj,jpk), sn_25h (jpi,jpj,jpk), sshn_25h(jpi,jpj)  ,     &
78         &      un_25h (jpi,jpj,jpk), vn_25h (jpi,jpj,jpk), wn_25h(jpi,jpj,jpk),     &
79         &      avt_25h(jpi,jpj,jpk), avm_25h(jpi,jpj,jpk),                      STAT=ierror )
80      IF( ierror > 0 ) THEN
81         CALL ctl_stop( 'dia_25h: unable to allocate ocean arrays' )   ;   RETURN
82      ENDIF
83      IF( ln_zdftke ) THEN             ! TKE physics
84         ALLOCATE( en_25h(jpi,jpj,jpk), STAT=ierror )
85         IF( ierror > 0 ) THEN
86            CALL ctl_stop( 'dia_25h: unable to allocate en_25h' )   ;   RETURN
87         ENDIF
88      ENDIF
89      IF( ln_zdfgls ) THEN             ! GLS physics
90         ALLOCATE( en_25h(jpi,jpj,jpk), rmxln_25h(jpi,jpj,jpk), STAT=ierror )
91         IF( ierror > 0 ) THEN
92            CALL ctl_stop( 'dia_25h: unable to allocate en_25h and rmxln_25h' )   ;   RETURN
93         ENDIF
94      ENDIF
95      ! ------------------------- !
96      ! 2 - Assign Initial Values !
97      ! ------------------------- !
98      cnt_25h = 1  ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible)
99      tn_25h  (:,:,:) = ts (:,:,:,jp_tem,Kbb)
100      sn_25h  (:,:,:) = ts (:,:,:,jp_sal,Kbb)
101      sshn_25h(:,:)   = ssh(:,:,Kbb)
102      un_25h  (:,:,:) = uu  (:,:,:,Kbb)
103      vn_25h  (:,:,:) = vv  (:,:,:,Kbb)
104      avt_25h (:,:,:) = avt (:,:,:)
105      avm_25h (:,:,:) = avm (:,:,:)
106      IF( ln_zdftke ) THEN
107         en_25h(:,:,:) = en(:,:,:)
108      ENDIF
109      IF( ln_zdfgls ) THEN
110         en_25h   (:,:,:) = en    (:,:,:)
111         rmxln_25h(:,:,:) = hmxl_n(:,:,:)
112      ENDIF
113#if defined key_si3
114      CALL ctl_stop('STOP', 'dia_25h not setup yet to do tidemean ice')
115#endif 
116      !
117   END SUBROUTINE dia_25h_init
118
119
120   SUBROUTINE dia_25h( kt, Kmm ) 
121      !!----------------------------------------------------------------------
122      !!                 ***  ROUTINE dia_25h  ***
123      !!         
124      !! ** Purpose :   Write diagnostics with M2/S2 tide removed
125      !!
126      !! ** Method  :   25hr mean outputs for shelf seas
127      !!----------------------------------------------------------------------
128      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
129      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
130      !!
131      INTEGER ::   ji, jj, jk
132      INTEGER                          ::   iyear0, nimonth0,iday0            ! start year,imonth,day
133      LOGICAL ::   ll_print = .FALSE.    ! =T print and flush numout
134      REAL(wp)                         ::   zsto, zout, zmax, zjulian, zmdi   ! local scalars
135      INTEGER                          ::   i_steps                           ! no of timesteps per hour
136      REAL(wp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                ! workspace
137      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d                              ! workspace
138      REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                             ! workspace
139      !!----------------------------------------------------------------------
140
141      ! 0. Initialisation
142      ! -----------------
143      ! Define frequency of summing to create 25 h mean
144      IF( MOD( 3600,NINT(rdt) ) == 0 ) THEN
145         i_steps = 3600/NINT(rdt)
146      ELSE
147         CALL ctl_stop('STOP', 'dia_wri_tide: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible')
148      ENDIF
149
150      ! local variable for debugging
151      ll_print = ll_print .AND. lwp
152
153      ! wn_25h could not be initialised in dia_25h_init, so we do it here instead
154      IF( kt == nn_it000 ) THEN
155         wn_25h(:,:,:) = ww(:,:,:)
156      ENDIF
157
158      ! Sum of 25 hourly instantaneous values to give a 25h mean from 24hours every day
159      IF( MOD( kt, i_steps ) == 0  .AND. kt /= nn_it000 ) THEN
160
161         IF (lwp) THEN
162              WRITE(numout,*) 'dia_wri_tide : Summing instantaneous hourly diagnostics at timestep ',kt
163              WRITE(numout,*) '~~~~~~~~~~~~ '
164         ENDIF
165
166         tn_25h  (:,:,:)     = tn_25h  (:,:,:) + ts (:,:,:,jp_tem,Kmm)
167         sn_25h  (:,:,:)     = sn_25h  (:,:,:) + ts (:,:,:,jp_sal,Kmm)
168         sshn_25h(:,:)       = sshn_25h(:,:)   + ssh(:,:,Kmm)
169         un_25h  (:,:,:)     = un_25h  (:,:,:) + uu  (:,:,:,Kmm)
170         vn_25h  (:,:,:)     = vn_25h  (:,:,:) + vv  (:,:,:,Kmm)
171         wn_25h  (:,:,:)     = wn_25h  (:,:,:) + ww  (:,:,:)
172         avt_25h (:,:,:)     = avt_25h (:,:,:) + avt (:,:,:)
173         avm_25h (:,:,:)     = avm_25h (:,:,:) + avm (:,:,:)
174         IF( ln_zdftke ) THEN
175            en_25h(:,:,:)    = en_25h  (:,:,:) + en(:,:,:)
176         ENDIF
177         IF( ln_zdfgls ) THEN
178            en_25h   (:,:,:) = en_25h   (:,:,:) + en    (:,:,:)
179            rmxln_25h(:,:,:) = rmxln_25h(:,:,:) + hmxl_n(:,:,:)
180         ENDIF
181         cnt_25h = cnt_25h + 1
182         !
183         IF (lwp) THEN
184            WRITE(numout,*) 'dia_tide : Summed the following number of hourly values so far',cnt_25h
185         ENDIF
186         !
187      ENDIF ! MOD( kt, i_steps ) == 0
188
189      ! Write data for 25 hour mean output streams
190      IF( cnt_25h == 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt /= nn_it000 ) THEN
191         !
192         IF(lwp) THEN
193            WRITE(numout,*) 'dia_wri_tide : Writing 25 hour mean tide diagnostics at timestep', kt
194            WRITE(numout,*) '~~~~~~~~~~~~ '
195         ENDIF
196         !
197         tn_25h  (:,:,:) = tn_25h  (:,:,:) * r1_25
198         sn_25h  (:,:,:) = sn_25h  (:,:,:) * r1_25
199         sshn_25h(:,:)   = sshn_25h(:,:)   * r1_25
200         un_25h  (:,:,:) = un_25h  (:,:,:) * r1_25
201         vn_25h  (:,:,:) = vn_25h  (:,:,:) * r1_25
202         wn_25h  (:,:,:) = wn_25h  (:,:,:) * r1_25
203         avt_25h (:,:,:) = avt_25h (:,:,:) * r1_25
204         avm_25h (:,:,:) = avm_25h (:,:,:) * r1_25
205         IF( ln_zdftke ) THEN
206            en_25h(:,:,:) = en_25h(:,:,:) * r1_25
207         ENDIF
208         IF( ln_zdfgls ) THEN
209            en_25h   (:,:,:) = en_25h   (:,:,:) * r1_25
210            rmxln_25h(:,:,:) = rmxln_25h(:,:,:) * r1_25
211         ENDIF
212         !
213         IF(lwp)  WRITE(numout,*) 'dia_wri_tide : Mean calculated by dividing 25 hour sums and writing output'
214         zmdi=1.e+20 !missing data indicator for masking
215         ! write tracers (instantaneous)
216         zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
217         CALL iom_put("temper25h", zw3d)   ! potential temperature
218         zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
219         CALL iom_put( "salin25h", zw3d  )   ! salinity
220         zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
221         IF( ll_wd ) THEN
222            CALL iom_put( "ssh25h", zw2d+ssh_ref )   ! sea surface
223         ELSE
224            CALL iom_put( "ssh25h", zw2d )   ! sea surface
225         ENDIF
226         ! Write velocities (instantaneous)
227         zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:))
228         CALL iom_put("vozocrtx25h", zw3d)    ! i-current
229         zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:))
230         CALL iom_put("vomecrty25h", zw3d  )   ! j-current
231         zw3d(:,:,:) = wn_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
232         CALL iom_put("vovecrtz25h", zw3d )   ! k-current
233         ! Write vertical physics
234         zw3d(:,:,:) = avt_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
235         CALL iom_put("avt25h", zw3d )   ! diffusivity
236         zw3d(:,:,:) = avm_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
237         CALL iom_put("avm25h", zw3d)   ! viscosity
238         IF( ln_zdftke ) THEN
239            zw3d(:,:,:) = en_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
240            CALL iom_put("tke25h", zw3d)   ! tke
241         ENDIF
242         IF( ln_zdfgls ) THEN
243            zw3d(:,:,:) = en_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
244            CALL iom_put("tke25h", zw3d)   ! tke
245            zw3d(:,:,:) = rmxln_25h(:,:,:)*wmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
246            CALL iom_put( "mxln25h",zw3d)
247         ENDIF
248         !
249         ! After the write reset the values to cnt=1 and sum values equal current value
250         tn_25h  (:,:,:) = ts (:,:,:,jp_tem,Kmm)
251         sn_25h  (:,:,:) = ts (:,:,:,jp_sal,Kmm)
252         sshn_25h(:,:)   = ssh(:,:,Kmm)
253         un_25h  (:,:,:) = uu  (:,:,:,Kmm)
254         vn_25h  (:,:,:) = vv  (:,:,:,Kmm)
255         wn_25h  (:,:,:) = ww  (:,:,:)
256         avt_25h (:,:,:) = avt (:,:,:)
257         avm_25h (:,:,:) = avm (:,:,:)
258         IF( ln_zdftke ) THEN
259            en_25h(:,:,:) = en(:,:,:)
260         ENDIF
261         IF( ln_zdfgls ) THEN
262            en_25h   (:,:,:) = en    (:,:,:)
263            rmxln_25h(:,:,:) = hmxl_n(:,:,:)
264         ENDIF
265         cnt_25h = 1
266         IF(lwp)  WRITE(numout,*) 'dia_wri_tide :   &
267            &    After 25hr mean write, reset sum to current value and cnt_25h to one for overlapping average', cnt_25h
268      ENDIF !  cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps * 24) == 0 .AND. kt .NE. nn_it000
269      !
270   END SUBROUTINE dia_25h 
271
272   !!======================================================================
273END MODULE dia25h
Note: See TracBrowser for help on using the repository browser.