source: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90 @ 8059

Last change on this file since 8059 was 8059, checked in by jgraham, 3 years ago

Merged branches required for AMM15 simulations, see ticket #1904.
Merged branches include:
branches/UKMO/CO6_KD490
branches/UKMO/CO6_Restartable_Tidal_Analysis
branches/UKMO/AMM15_v3_6_STABLE

File size: 14.2 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 diainsitutem, ONLY: rinsitu_t, theta2t
11   USE in_out_manager  ! I/O units
12   USE iom             ! I/0 library
13   USE wrk_nemo        ! working arrays
14#if defined key_zdftke 
15   USE zdftke, ONLY: en
16#endif
17   USE zdf_oce, ONLY: avt, avm
18#if defined key_zdfgls
19   USE zdfgls, ONLY: mxln
20   USE zdf_oce, ONLY: en
21#endif
22
23   IMPLICIT NONE
24   PRIVATE
25
26   LOGICAL , PUBLIC ::   ln_dia25h     !:  25h mean output
27   PUBLIC   dia_25h_init               ! routine called by nemogcm.F90
28   PUBLIC   dia_25h                    ! routine called by diawri.F90
29
30  !! * variables for calculating 25-hourly means
31   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_25h  , sn_25h, rinsitu_t_25h 
32   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_25h 
33   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   un_25h  , vn_25h  , wn_25h
34   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   avt_25h , avm_25h
35#if defined key_zdfgls || key_zdftke
36   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   en_25h
37#endif
38#if defined key_zdfgls 
39   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   rmxln_25h
40#endif
41   INTEGER, SAVE :: cnt_25h     ! Counter for 25 hour means
42
43
44
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   SUBROUTINE dia_25h_init 
53      !!---------------------------------------------------------------------------
54      !!                  ***  ROUTINE dia_25h_init  ***
55      !!     
56      !! ** Purpose: Initialization of 25h mean namelist
57      !!       
58      !! ** Method : Read namelist
59      !!   History
60      !!   3.6  !  08-14  (E. O'Dea) Routine to initialize dia_25h
61      !!---------------------------------------------------------------------------
62      !!
63      INTEGER ::   ios                 ! Local integer output status for namelist read
64      INTEGER ::   ierror              ! Local integer for memory allocation
65      !
66      NAMELIST/nam_dia25h/ ln_dia25h
67      !!----------------------------------------------------------------------
68      !
69      REWIND ( numnam_ref )              ! Read Namelist nam_dia25h in reference namelist : 25hour mean diagnostics
70      READ   ( numnam_ref, nam_dia25h, IOSTAT=ios, ERR= 901 )
71901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_dia25h in reference namelist', lwp )
72
73      REWIND( numnam_cfg )              ! Namelist nam_dia25h in configuration namelist  25hour diagnostics
74      READ  ( numnam_cfg, nam_dia25h, IOSTAT = ios, ERR = 902 )
75902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_dia25h in configuration namelist', lwp )
76      IF(lwm) WRITE ( numond, nam_dia25h )
77
78      IF(lwp) THEN                   ! Control print
79         WRITE(numout,*)
80         WRITE(numout,*) 'dia_25h_init : Output 25 hour mean diagnostics'
81         WRITE(numout,*) '~~~~~~~~~~~~'
82         WRITE(numout,*) 'Namelist nam_dia25h : set 25h outputs '
83         WRITE(numout,*) 'Switch for 25h diagnostics (T) or not (F)  ln_dia25h  = ', ln_dia25h
84      ENDIF
85      IF( .NOT. ln_dia25h )   RETURN
86      ! ------------------- !
87      ! 1 - Allocate memory !
88      ! ------------------- !
89      ALLOCATE( tn_25h(jpi,jpj,jpk), STAT=ierror )
90      IF( ierror > 0 ) THEN
91         CALL ctl_stop( 'dia_25h: unable to allocate tn_25h' )   ;   RETURN
92      ENDIF
93      ALLOCATE( sn_25h(jpi,jpj,jpk), STAT=ierror )
94      IF( ierror > 0 ) THEN
95         CALL ctl_stop( 'dia_25h: unable to allocate sn_25h' )   ;   RETURN
96      ENDIF
97      ALLOCATE( rinsitu_t_25h(jpi,jpj,jpk), STAT=ierror )
98      IF( ierror > 0 ) THEN
99         CALL ctl_stop( 'dia_25h: unable to allocate rinsitu_t_25h' )   ;   RETURN
100      ENDIF
101      ALLOCATE( un_25h(jpi,jpj,jpk), STAT=ierror )
102      IF( ierror > 0 ) THEN
103         CALL ctl_stop( 'dia_25h: unable to allocate un_25h' )   ;   RETURN
104      ENDIF
105      ALLOCATE( vn_25h(jpi,jpj,jpk), STAT=ierror )
106      IF( ierror > 0 ) THEN
107         CALL ctl_stop( 'dia_25h: unable to allocate vn_25h' )   ;   RETURN
108      ENDIF
109      ALLOCATE( wn_25h(jpi,jpj,jpk), STAT=ierror )
110      IF( ierror > 0 ) THEN
111         CALL ctl_stop( 'dia_25h: unable to allocate wn_25h' )   ;   RETURN
112      ENDIF
113      ALLOCATE( avt_25h(jpi,jpj,jpk), STAT=ierror )
114      IF( ierror > 0 ) THEN
115         CALL ctl_stop( 'dia_25h: unable to allocate avt_25h' )   ;   RETURN
116      ENDIF
117      ALLOCATE( avm_25h(jpi,jpj,jpk), STAT=ierror )
118      IF( ierror > 0 ) THEN
119         CALL ctl_stop( 'dia_25h: unable to allocate avm_25h' )   ;   RETURN
120      ENDIF
121# if defined key_zdfgls || defined key_zdftke
122      ALLOCATE( en_25h(jpi,jpj,jpk), STAT=ierror )
123      IF( ierror > 0 ) THEN
124         CALL ctl_stop( 'dia_25h: unable to allocate en_25h' )   ;   RETURN
125      ENDIF
126#endif
127# if defined key_zdfgls 
128      ALLOCATE( rmxln_25h(jpi,jpj,jpk), STAT=ierror )
129      IF( ierror > 0 ) THEN
130         CALL ctl_stop( 'dia_25h: unable to allocate rmxln_25h' )   ;   RETURN
131      ENDIF
132#endif
133      ALLOCATE( sshn_25h(jpi,jpj), STAT=ierror )
134      IF( ierror > 0 ) THEN
135         CALL ctl_stop( 'dia_25h: unable to allocate sshn_25h' )   ;   RETURN
136      ENDIF
137      ! ------------------------- !
138      ! 2 - Assign Initial Values !
139      ! ------------------------- !
140      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)
141      tn_25h(:,:,:) = tsb(:,:,:,jp_tem)
142      sn_25h(:,:,:) = tsb(:,:,:,jp_sal)
143      CALL theta2t
144      rinsitu_t_25h(:,:,:) = rinsitu_t(:,:,:)
145      sshn_25h(:,:) = sshb(:,:)
146      un_25h(:,:,:) = ub(:,:,:)
147      vn_25h(:,:,:) = vb(:,:,:)
148      wn_25h(:,:,:) = wn(:,:,:)
149      avt_25h(:,:,:) = avt(:,:,:)
150      avm_25h(:,:,:) = avm(:,:,:)
151# if defined key_zdfgls || defined key_zdftke
152         en_25h(:,:,:) = en(:,:,:)
153#endif
154# if defined key_zdfgls
155         rmxln_25h(:,:,:) = mxln(:,:,:)
156#endif
157#if defined key_lim3 || defined key_lim2
158         CALL ctl_stop('STOP', 'dia_25h not setup yet to do tidemean ice')
159#endif 
160
161      ! -------------------------- !
162      ! 3 - Return to dia_wri      !
163      ! -------------------------- !
164
165
166   END SUBROUTINE dia_25h_init
167
168
169   SUBROUTINE dia_25h( kt ) 
170      !!----------------------------------------------------------------------
171      !!                 ***  ROUTINE dia_25h  ***
172      !!         
173      !!
174      !!--------------------------------------------------------------------
175      !!                   
176      !! ** Purpose :   Write diagnostics with M2/S2 tide removed
177      !!
178      !! ** Method  :   
179      !!      25hr mean outputs for shelf seas
180      !!
181      !! History :
182      !!   ?.0  !  07-04  (A. Hines) New routine, developed from dia_wri_foam
183      !!   3.4  !  02-13  (J. Siddorn) Routine taken from old dia_wri_foam
184      !!   3.6  !  08-14  (E. O'Dea) adapted for VN3.6
185      !!----------------------------------------------------------------------
186      !! * Modules used
187
188      IMPLICIT NONE
189
190      !! * Arguments
191      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
192
193
194      !! * Local declarations
195      INTEGER ::   ji, jj, jk
196
197      LOGICAL ::   ll_print = .FALSE.    ! =T print and flush numout
198      REAL(wp)                         ::   zsto, zout, zmax, zjulian, zdt, zmdi  ! temporary reals
199      INTEGER                          ::   i_steps                               ! no of timesteps per hour
200      REAL(wp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                    ! temporary workspace
201      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d                                  ! temporary workspace
202      REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                                 ! temporary workspace
203      INTEGER                          ::   iyear0, nimonth0,iday0                ! start year,imonth,day
204
205      !!----------------------------------------------------------------------
206
207      ! 0. Initialisation
208      ! -----------------
209      ! Define frequency of summing to create 25 h mean
210      zdt = rdt
211      IF( nacc == 1 ) zdt = rdtmin
212
213      IF( MOD( 3600,INT(zdt) ) == 0 ) THEN
214         i_steps = 3600/INT(zdt)
215      ELSE
216         CALL ctl_stop('STOP', 'dia_wri_tide: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible')
217      ENDIF
218
219#if defined key_lim3 || defined key_lim2
220      CALL ctl_stop('STOP', 'dia_wri_tide not setup yet to do tidemean ice')
221#endif
222
223      ! local variable for debugging
224      ll_print = ll_print .AND. lwp
225
226      ! Sum of 25 hourly instantaneous values to give a 25h mean from 24hours
227      ! every day
228      IF( MOD( kt, i_steps ) == 0  .and. kt .ne. nn_it000 ) THEN
229
230         IF (lwp) THEN
231              WRITE(numout,*) 'dia_wri_tide : Summing instantaneous hourly diagnostics at timestep ',kt
232              WRITE(numout,*) '~~~~~~~~~~~~ '
233         ENDIF
234
235         tn_25h(:,:,:)        = tn_25h(:,:,:) + tsn(:,:,:,jp_tem)
236         sn_25h(:,:,:)        = sn_25h(:,:,:) + tsn(:,:,:,jp_sal)
237         CALL theta2t
238         rinsitu_t_25h(:,:,:)  = rinsitu_t_25h(:,:,:) + rinsitu_t(:,:,:)
239         sshn_25h(:,:)        = sshn_25h(:,:) + sshn (:,:)
240         un_25h(:,:,:)        = un_25h(:,:,:) + un(:,:,:)
241         vn_25h(:,:,:)        = vn_25h(:,:,:) + vn(:,:,:)
242         wn_25h(:,:,:)        = wn_25h(:,:,:) + wn(:,:,:)
243         avt_25h(:,:,:)       = avt_25h(:,:,:) + avt(:,:,:)
244         avm_25h(:,:,:)       = avm_25h(:,:,:) + avm(:,:,:)
245# if defined key_zdfgls || defined key_zdftke
246         en_25h(:,:,:)        = en_25h(:,:,:) + en(:,:,:)
247#endif
248# if defined key_zdfgls
249         rmxln_25h(:,:,:)      = rmxln_25h(:,:,:) + mxln(:,:,:)
250#endif
251         cnt_25h = cnt_25h + 1
252
253         IF (lwp) THEN
254            WRITE(numout,*) 'dia_tide : Summed the following number of hourly values so far',cnt_25h
255         ENDIF
256
257      ENDIF ! MOD( kt, i_steps ) == 0
258
259         ! Write data for 25 hour mean output streams
260      IF( cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN
261
262            IF(lwp) THEN
263               WRITE(numout,*) 'dia_wri_tide : Writing 25 hour mean tide diagnostics at timestep', kt
264               WRITE(numout,*) '~~~~~~~~~~~~ '
265            ENDIF
266
267            tn_25h(:,:,:)        = tn_25h(:,:,:) / 25.0_wp
268            sn_25h(:,:,:)        = sn_25h(:,:,:) / 25.0_wp
269            rinsitu_t_25h(:,:,:)  = rinsitu_t_25h(:,:,:) / 25.0_wp
270            sshn_25h(:,:)        = sshn_25h(:,:) / 25.0_wp
271            un_25h(:,:,:)        = un_25h(:,:,:) / 25.0_wp
272            vn_25h(:,:,:)        = vn_25h(:,:,:) / 25.0_wp
273            wn_25h(:,:,:)        = wn_25h(:,:,:) / 25.0_wp
274            avt_25h(:,:,:)       = avt_25h(:,:,:) / 25.0_wp
275            avm_25h(:,:,:)       = avm_25h(:,:,:) / 25.0_wp
276# if defined key_zdfgls || defined key_zdftke
277            en_25h(:,:,:)        = en_25h(:,:,:) / 25.0_wp
278#endif
279# if defined key_zdfgls
280            rmxln_25h(:,:,:)       = rmxln_25h(:,:,:) / 25.0_wp
281#endif
282
283            IF (lwp)  WRITE(numout,*) 'dia_wri_tide : Mean calculated by dividing 25 hour sums and writing output'
284            zmdi=1.e+20 !missing data indicator for masking
285            ! write tracers (instantaneous)
286            zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
287            CALL iom_put("temper25h", zw3d)   ! potential temperature
288            CALL theta2t                                                                    ! calculate insitu temp
289            zw3d(:,:,:) = rinsitu_t_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
290            CALL iom_put("tempis25h", zw3d)   ! in-situ temperature
291            zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
292            CALL iom_put( "salin25h", zw3d  )   ! salinity
293            zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
294            CALL iom_put( "ssh25h", zw2d )   ! sea surface
295
296
297            ! Write velocities (instantaneous)
298            zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:))
299            CALL iom_put("vozocrtx25h", zw3d)    ! i-current
300            zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:))
301            CALL iom_put("vomecrty25h", zw3d  )   ! j-current
302
303            zw3d(:,:,:) = wn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
304            CALL iom_put("vomecrtz25h", zw3d )   ! k-current
305            zw3d(:,:,:) = avt_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
306            CALL iom_put("avt25h", zw3d )   ! diffusivity
307            zw3d(:,:,:) = avm_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
308            CALL iom_put("avm25h", zw3d)   ! viscosity
309#if defined key_zdftke || defined key_zdfgls 
310!           zw3d(:,:,:) = en_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
311!           CALL iom_put("tke25h", zw3d)   ! tke
312#endif
313#if defined key_zdfgls 
314!           zw3d(:,:,:) = rmxln_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
315!           CALL iom_put( "mxln25h",zw3d)
316#endif
317
318            ! After the write reset the values to cnt=1 and sum values equal current value
319            tn_25h(:,:,:) = tsn(:,:,:,jp_tem)
320            sn_25h(:,:,:) = tsn(:,:,:,jp_sal)
321            CALL theta2t
322            rinsitu_t_25h(:,:,:) = rinsitu_t(:,:,:)
323            sshn_25h(:,:) = sshn (:,:)
324            un_25h(:,:,:) = un(:,:,:)
325            vn_25h(:,:,:) = vn(:,:,:)
326            wn_25h(:,:,:) = wn(:,:,:)
327            avt_25h(:,:,:) = avt(:,:,:)
328            avm_25h(:,:,:) = avm(:,:,:)
329# if defined key_zdfgls || defined key_zdftke
330            en_25h(:,:,:) = en(:,:,:)
331#endif
332# if defined key_zdfgls
333            rmxln_25h(:,:,:) = mxln(:,:,:)
334#endif
335            cnt_25h = 1
336            IF (lwp)  WRITE(numout,*) 'dia_wri_tide : After 25hr mean write, reset sum to current value and cnt_25h to one for overlapping average',cnt_25h
337
338      ENDIF !  cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps * 24) == 0 .AND. kt .NE. nn_it000
339
340
341   END SUBROUTINE dia_25h 
342
343   !!======================================================================
344END MODULE dia25h
Note: See TracBrowser for help on using the repository browser.