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.
diaopfoam.F90 in branches/UKMO/AMM15_v3_6_STABLE_package_collate_PS44/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_PS44/NEMOGCM/NEMO/OPA_SRC/DIA/diaopfoam.F90 @ 12536

Last change on this file since 12536 was 12536, checked in by petesykes, 4 years ago

added diaopfoam T0 diagnostics
File size: 18.7 KB
Line 
1MODULE diaopfoam
2   !!======================================================================
3   !!                       ***  MODULE  diaopfoam  ***
4   !! Output stream for operational use
5   !!======================================================================
6   !! History :  3.6  !  2016  (P Sykes)  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   USE diatmb
15   USE diurnal_bulk
16   USE cool_skin
17   USE ioipsl
18
19
20   IMPLICIT NONE
21   PRIVATE
22
23   LOGICAL , PUBLIC ::   ln_diaopfoam     !: Diaopfoam output
24   LOGICAL , PUBLIC ::   ln_diaopfoam_Tzero     !: Diaopfoam first time step output
25   PUBLIC   dia_diaopfoam_init            ! routine called by nemogcm.F90
26   PUBLIC   dia_diaopfoam                 ! routine called by diawri.F90
27   PUBLIC   calc_max_cur                  ! routine called by diaopfoam.F90
28
29   !! * Substitutions
30#  include "domzgr_substitute.h90"
31
32   !!----------------------------------------------------------------------
33   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
34   !! $Id$
35   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
36   !!----------------------------------------------------------------------
37CONTAINS
38
39   SUBROUTINE dia_diaopfoam_init
40      !!---------------------------------------------------------------------------
41      !!                  ***  ROUTINE dia_wri_diaop_init  ***
42      !!     
43      !! ** Purpose: Initialization of diaopfoam namelist
44      !!       
45      !! ** Method : Read namelist
46      !!   History
47      !!   3.4  !  03-14  (P. Sykes) Routine to initialize dia_wri_diaop
48      !!---------------------------------------------------------------------------
49      !!
50      INTEGER            ::   ios      ! Local integer output status for namelist read
51      INTEGER            ::   ierror   ! local integer
52      !!
53      NAMELIST/nam_diadiaop/ ln_diaopfoam,ln_diaopfoam_Tzero
54      !!
55      !!----------------------------------------------------------------------
56      !
57      ln_diaopfoam = .false.         ! default value for diaopfoam stream
58      ln_diaopfoam_Tzero = .false.         ! default value for diaopfoam Tzero stream
59      REWIND ( numnam_ref )              ! Read Namelist nam_diadiaop in reference namelist : 3D hourly diagnostics
60      READ   ( numnam_ref, nam_diadiaop, IOSTAT=ios, ERR= 901 )
61901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadiaop in reference namelist', lwp )
62
63      REWIND( numnam_cfg )              ! Namelist nam_diadiaop in configuration namelist  3D hourly diagnostics
64      READ  ( numnam_cfg, nam_diadiaop, IOSTAT = ios, ERR = 902 )
65902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diadiaop in configuration namelist', lwp )
66      IF(lwm) WRITE ( numond, nam_diadiaop )
67      !
68      IF(lwp) THEN                   ! Control print
69         WRITE(numout,*)
70         WRITE(numout,*) 'dia_diaopfoam_init : Output Diaopfoam Diagnostics'
71         WRITE(numout,*) '~~~~~~~~~~~~'
72         WRITE(numout,*) '   Namelist nam_diadiaop : set diaopfoam outputs '
73         WRITE(numout,*) '      Switch for diaopfoam diagnostics (T) or not (F)                 ln_diaopfoam        = ', ln_diaopfoam
74         WRITE(numout,*) '      Switch for diaopfoam first timestep diagnostics (T) or not (F)  ln_diaopfoam_Tzero  = ', ln_diaopfoam_Tzero
75      ENDIF
76    END SUBROUTINE dia_diaopfoam_init
77
78    SUBROUTINE dia_diaopfoam( kt )
79      !!----------------------------------------------------------------------
80      !!                 ***  ROUTINE dia_diaopfoam  ***
81      !! ** Purpose :   Write 3D hourly diagnostics for operational use
82      !!
83      !!
84      !! History :
85      !!   3.6  !  11-16  (P. Sykes)
86      !!         
87      !!--------------------------------------------------------------------
88      IMPLICIT NONE
89
90      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
91
92      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
93      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
94      REAL(wp) :: zmdi
95      REAL(wp), POINTER, DIMENSION(:,:)   :: zwu
96      REAL(wp), POINTER, DIMENSION(:,:)   :: zwv
97      REAL(wp), POINTER, DIMENSION(:,:)   :: zwz
98
99      CALL wrk_alloc( jpi , jpj      , zwu )
100      CALL wrk_alloc( jpi , jpj      , zwv )
101      CALL wrk_alloc( jpi , jpj      , zwz )
102
103      zmdi=1.e+20                               !  missing data indicator for masking
104
105      ! Diaopfoam stream if needed
106      IF (ln_diaopfoam) THEN
107         IF ( kt .eq. nn_it000 .AND. ln_diaopfoam_Tzero ) THEN
108            IF(lwp) WRITE(numout,*) 'diaopfoam: writing T0 at kt = ', kt
109            CALL dia_diaopfoam_zero( 'Tzero', kt )
110         ENDIF
111         
112         CALL theta2t
113         CALL iom_put( "insitut_op" , rinsitu_t(:,:,:)                      )    ! insitu temperature
114         CALL iom_put( "toce_op"   , tsn(:,:,:,jp_tem)                     )    ! temperature
115         CALL iom_put( "soce_op"   , tsn(:,:,:,jp_sal)                     )    ! salinity
116         IF (ln_diurnal) THEN
117            CALL iom_put( "sst_wl_op"   , x_dsst                             )    ! warm layer
118            CALL iom_put( "sst_cs_op"   , x_csdsst                           )    ! cool skin
119         ENDIF
120         zw2d(:,:)=sshn(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
121         CALL iom_put( "ssh_op"  , zw2d(:,:)                          ) ! sea surface height
122         CALL iom_put( "uoce_op"   , un                                    )    ! i-current     
123         CALL iom_put( "voce_op"   , vn                                    )    ! j-current
124         !CALL iom_put( "woce_op"   , wn                                    )    ! k-current
125         CALL calc_max_cur(zwu,zwv,zwz,zmdi)
126         CALL iom_put( "maxu" , zwu                                     ) ! max u current
127         CALL iom_put( "maxv" , zwv                                     ) ! max v current
128         CALL iom_put( "maxz" , zwz                                     ) ! max current depth
129      ENDIF
130   END SUBROUTINE dia_diaopfoam
131
132   SUBROUTINE calc_max_cur(zmax_u, zmax_v, zmax_z, inmdi)
133      !!---------------------------------------------------------------------
134      !!                    ***  ROUTINE calc_max_cur  ***
135      !!
136      !! ** Purpose :   To locate within the water column the magnitude and
137      !!                vertical location of the strongest horizontal
138      !!                current. The vertical component is ignored since it
139      !!                is an order of magnitude smaller than the horizontal
140      !!                flow, in general.
141      !!
142      !! ** Method :    A. Map U,V to T-grid.
143      !!                B. Calculate the magnitude of the current for every
144      !!                   grid cell.
145      !!                C. Locate the vertical index using FORTRAN's builtin
146      !!                   MAXLOC function.
147      !!                D. Copy the U,V,Z components of the relevant
148      !!                   indices.
149      !!
150      !! ** Returns :   Value of u, v component and depth of maximum
151      !!                horizontal current on T-grid.
152      !!
153      !!---------------------------------------------------------------------
154      IMPLICIT NONE
155      REAL(wp), DIMENSION(jpi, jpj), INTENT(  OUT) :: zmax_u, zmax_v, zmax_z
156      REAL(wp), DIMENSION(jpi, jpj, jpk)           :: zmax_u_t, zmax_v_t
157      REAL(wp), DIMENSION(jpi, jpj, jpk)           :: zcmag
158      REAL(wp), DIMENSION(jpi, jpj)                :: zmaxk
159      REAL(wp),                      INTENT(IN   ) :: inmdi
160      INTEGER :: ji, jj, jk
161
162      ! Initialise output arrays
163      zmax_u(:, :) = inmdi
164      zmax_v(:, :) = inmdi
165      zmax_z(:, :) = inmdi
166      ! Map to T-grid
167      zmax_u_t(:, :, :) = 0._wp
168      zmax_v_t(:, :, :) = 0._wp
169      DO jk = 1,jpk
170         DO jj = 2,jpj
171            DO ji = 2,jpi
172               zmax_u_t(ji, jj, jk) = 0.5 * (un(ji, jj, jk) + un(ji-1, jj, jk))
173               zmax_v_t(ji, jj, jk) = 0.5 * (vn(ji, jj, jk) + vn(ji, jj-1, jk))
174            END DO
175         END DO
176      END DO
177      ! Calculate absolute velocity
178      zcmag = sqrt((zmax_u_t)**2 + (zmax_v_t)**2)
179      ! Find max. current
180      zmaxk = maxloc(zcmag, dim=3)
181      ! Output values
182      DO jj = 1,jpj
183         DO ji = 1,jpi
184             zmax_u(ji, jj) = zmax_u_t(ji, jj, INT(zmaxk(ji, jj)))
185             zmax_v(ji, jj) = zmax_v_t(ji, jj, INT(zmaxk(ji, jj)))
186             zmax_z(ji, jj) = fsdept(ji, jj, INT(zmaxk(ji, jj)))
187         END DO
188      END DO     
189   END SUBROUTINE calc_max_cur
190
191   SUBROUTINE dia_diaopfoam_zero( cdfile_name, kt )
192      !!---------------------------------------------------------------------
193      !!                 ***  ROUTINE dia_diaopfoam_zero  ***
194      !!       
195      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
196      !!      the instantaneous ocean state at the first tiome step.
197      !!
198      !! ** Method  :   NetCDF files using ioipsl
199      !!----------------------------------------------------------------------
200      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
201      INTEGER           , INTENT( in ) ::   kt               ! ocean time-step index
202      !!
203      CHARACTER (len=32) :: clhstnam
204      CHARACTER (len=40) :: clop
205      INTEGER  ::   iimi, iima, ipk, ijmi, ijma          ! local integers
206      INTEGER  ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
207      INTEGER  ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
208      INTEGER  ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
209      INTEGER  ::   id_i , nz_i, nh_i       
210      INTEGER, DIMENSION(1) ::   idex                    ! local workspace
211      INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
212      INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
213      INTEGER :: ierr
214      INTEGER :: jkbot, jj, ji
215      REAL(wp) :: zsto, zout, zmax, zjulian, zdt
216      REAL(wp) :: zmdi
217      REAL(wp), POINTER, DIMENSION(:,:)   :: zwu
218      REAL(wp), POINTER, DIMENSION(:,:)   :: zwv
219      REAL(wp), POINTER, DIMENSION(:,:)   :: zwz
220      REAL(wp), DIMENSION(jpi,jpj) :: zw2d           ! 2D workspace
221      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
222      REAL(wp), DIMENSION(jpi,jpj) :: z2d
223
224      !!----------------------------------------------------------------------
225
226      ! -----------------------------------------------------------------
227      ! 0. Allocations
228      ! -----------------------------------------------------------------
229      ierr = 0
230      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
231         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
232         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr )
233      IF( lk_mpp )   CALL mpp_sum( ierr )
234         IF( ierr /= 0 ) THEN
235            CALL ctl_stop('dia_diaopfoam_zero: failed to allocate arrays')
236            RETURN
237      ENDIF
238      CALL wrk_alloc( jpi , jpj      , zwu )
239      CALL wrk_alloc( jpi , jpj      , zwv )
240      CALL wrk_alloc( jpi , jpj      , zwz )
241
242      zmdi=1.e+20                               !  missing data indicator for masking
243
244      ! -----------------------------------------------------------------
245      ! 1. Define NETCDF files and fields at beginning of first time step
246      ! -----------------------------------------------------------------
247
248      ! Define indices of the horizontal output zoom and vertical limit storage
249      iimi = 1      ;      iima = jpi
250      ijmi = 1      ;      ijma = jpj
251      ipk = jpk
252
253      ! Define frequency of output and means
254      zdt  = rdt
255      zsto = rdt
256      clop = "inst(x)"           ! no use of the mask value (require less cpu time)
257      zout = rdt
258      zmax = ( nitend - nit000 + 1 ) * zdt
259
260      ! Compute julian date from starting date of the run
261      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )         ! time axis
262      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
263
264      ! Define the T grid FILE ( nid_T )
265      clhstnam = TRIM(cdfile_name)//".grid_T"
266      IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
267      CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
268         &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
269         &          nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
270      CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
271         &           "m", ipk, gdept_1d, nz_T, "down" )
272      !                                                            ! Index of ocean points
273      CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
274      CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
275
276      ! Define the U grid FILE ( nid_U )
277      clhstnam = TRIM(cdfile_name)//".grid_U"
278      IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
279      CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
280         &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
281         &          nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
282      CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdepu
283         &           "m", ipk, gdept_1d, nz_U, "down" )
284      !                                                            ! Index of ocean points
285      CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
286      CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
287
288      ! Define the V grid FILE ( nid_V )
289      clhstnam = TRIM(cdfile_name)//".grid_V"
290      IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
291      CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
292         &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
293         &          nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
294      CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdepv
295         &          "m", ipk, gdept_1d, nz_V, "down" )
296      !                                                            ! Index of ocean points
297      CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
298      CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
299
300
301      ! -----------------------------------------------------------------
302      ! 2. Declare all the output fields as NETCDF variables
303      ! -----------------------------------------------------------------
304
305      !                                                                                     !!! nid_T : 3D
306      !CALL histdef( nid_T, "votempis", "Insitu Temperature" , "C"      ,       &  !
307      !   &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
308      CALL histdef( nid_T, "votemper", "Temperature"         , "C"     ,       &  ! tn
309         &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
310      CALL histdef( nid_T, "vosaline", "Salinity"            , "PSU"   ,       & ! sn
311         &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
312      CALL histdef( nid_T, "sossheig", "Sea Surface Height"  , "m"     ,       &  ! sshn
313         &          jpi, jpj, nh_T, 1  , 1, 1  , nz_T, 32, clop, zsto, zout )
314      CALL histdef( nid_T, "votempis", "Insitu Temperature"  , "C"     ,       &  ! rinsitu_t
315         &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
316      CALL histdef( nid_T, "maxu"    , "Max Zonal Current"   , "m/s"     ,       &  ! zwu
317         &          jpi, jpj, nh_T, 1  , 1, 1  , nz_T, 32, clop, zsto, zout )
318      CALL histdef( nid_T, "maxv"    , "Max Meridional Current" , "m/s"     ,    &  ! zwv
319         &          jpi, jpj, nh_T, 1  , 1, 1  , nz_T, 32, clop, zsto, zout )
320      CALL histdef( nid_T, "maxz"    , "Max Current Depth"   , "m/s"     ,       &  ! zwz
321         &          jpi, jpj, nh_T, 1  , 1, 1  , nz_T, 32, clop, zsto, zout )
322      CALL histdef( nid_T, "sbt"       , "Bottom Temperature"  , "C"     ,       &  ! sbt
323         &          jpi, jpj, nh_T, 1  , 1, 1  , nz_T, 32, clop, zsto, zout )
324      CALL histend( nid_T, snc4chunks=snc4set )
325
326      !                                                                                     !!! nid_U : 3D
327      CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
328         &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
329      CALL histend( nid_U, snc4chunks=snc4set )
330
331      !                                                                                     !!! nid_V : 3D
332      CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
333         &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
334      CALL histend( nid_V, snc4chunks=snc4set )
335
336
337      ! -----------------------------------------------------------------
338      ! 3. Write the data
339      ! -----------------------------------------------------------------
340
341      idex(1) = 1
342      CALL histwrite( nid_T, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature
343      CALL histwrite( nid_T, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity
344      CALL histwrite( nid_T, "sossheig", kt, sshn             , jpi*jpj    , idex )    ! sea surface height
345      CALL theta2t
346      CALL histwrite( nid_T, "votempis", kt, rinsitu_t(:,:,:) , jpi*jpj*jpk, idex )    ! now insitu-temperature
347      CALL calc_max_cur(zwu,zwv,zwz,zmdi)
348      CALL lbc_lnk( zwu, 'T', 1. )
349      CALL lbc_lnk( zwv, 'T', 1. )
350      CALL lbc_lnk( zwz, 'T', 1. )
351      CALL histwrite( nid_T, "maxu"    , kt, zwu              , jpi*jpj    , idex )    ! max u-current
352      CALL histwrite( nid_T, "maxv"    , kt, zwv              , jpi*jpj    , idex )    ! max v-current
353      CALL histwrite( nid_T, "maxz"    , kt, zwz              , jpi*jpj    , idex )    ! max current depth
354      DO jj = 1, jpj
355         DO ji = 1, jpi
356            jkbot = mbkt(ji,jj)
357            z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem)
358         END DO
359      END DO
360      CALL histwrite( nid_T, "sbt"     , kt, z2d              , jpi*jpj    , idex )    ! sbt
361      ! U file
362      CALL histwrite( nid_U, "vozocrtx", kt, un              , jpi*jpj*jpk, idex )     ! now i-velocity
363      ! V file
364      CALL histwrite( nid_V, "vomecrty", kt, vn              , jpi*jpj*jpk, idex )     ! now j-velocity
365
366
367      ! -----------------------------------------------------------------
368      ! 4. Close the files
369      ! -----------------------------------------------------------------
370
371      CALL histclo( nid_T )
372      CALL histclo( nid_U )
373      CALL histclo( nid_V )
374
375   END SUBROUTINE dia_diaopfoam_zero
376
377END MODULE diaopfoam
Note: See TracBrowser for help on using the repository browser.