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 branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

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

Last change on this file since 8561 was 8561, checked in by jgraham, 7 years ago

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