source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/STATION_ASF/MY_SRC/diawri.F90 @ 13159

Last change on this file since 13159 was 12489, checked in by davestorkey, 8 months ago

Preparation for new timestepping scheme #2390.
Main changes:

  1. Initial euler timestep now handled in stp and not in TRA/DYN routines.
  2. Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
  3. Renaming of a few miscellaneous parameters, eg. atfp → rn_atfp (namelist parameter used everywhere) and rau0 → rho0.

This version gives bit-comparable results to the previous version of the trunk.

File size: 19.3 KB
Line 
1MODULE diawri
2   !!======================================================================
3   !!                     ***  MODULE  diawri  ***
4   !! Ocean diagnostics :  write ocean output files
5   !!=====================================================================
6   !! History :  OPA  ! 1991-03  (M.-A. Foujols)  Original code
7   !!            4.0  ! 1991-11  (G. Madec)
8   !!                 ! 1992-06  (M. Imbard)  correction restart file
9   !!                 ! 1992-07  (M. Imbard)  split into diawri and rstwri
10   !!                 ! 1993-03  (M. Imbard)  suppress writibm
11   !!                 ! 1998-01  (C. Levy)  NETCDF format using ioipsl INTERFACE
12   !!                 ! 1999-02  (E. Guilyardi)  name of netCDF files + variables
13   !!            8.2  ! 2000-06  (M. Imbard)  Original code (diabort.F)
14   !!   NEMO     1.0  ! 2002-06  (A.Bozec, E. Durand)  Original code (diainit.F)
15   !!             -   ! 2002-09  (G. Madec)  F90: Free form and module
16   !!             -   ! 2002-12  (G. Madec)  merge of diabort and diainit, F90
17   !!                 ! 2005-11  (V. Garnier) Surface pressure gradient organization
18   !!            3.2  ! 2008-11  (B. Lemaire) creation from old diawri
19   !!            3.7  ! 2014-01  (G. Madec) remove eddy induced velocity from no-IOM output
20   !!                 !                     change name of output variables in dia_wri_state
21   !!----------------------------------------------------------------------
22
23   !!----------------------------------------------------------------------
24   !!   dia_wri       : create the standart output files
25   !!   dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
26   !!----------------------------------------------------------------------
27   USE oce            ! ocean dynamics and tracers
28   USE dom_oce        ! ocean space and time domain
29   USE phycst         ! physical constants
30   USE dianam         ! build name of file (routine)
31   USE sbc_oce        ! Surface boundary condition: ocean fields
32   !
33   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
34   USE in_out_manager ! I/O manager
35   USE iom            !
36   USE ioipsl         !
37#if defined key_si3
38   USE ice
39   USE icewri
40#endif
41   USE lib_mpp         ! MPP library
42   USE timing          ! preformance summary
43
44   IMPLICIT NONE
45   PRIVATE
46
47   PUBLIC   dia_wri                 ! routines called by step.F90
48   PUBLIC   dia_wri_state
49   PUBLIC   dia_wri_alloc           ! Called by nemogcm module
50
51   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
52   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
53   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
54   INTEGER ::   ndex(1)                              ! ???
55   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
56
57   !!----------------------------------------------------------------------
58   !! NEMO/SAS 4.0 , NEMO Consortium (2018)
59   !! $Id: diawri.F90 10425 2018-12-19 21:54:16Z smasson $
60   !! Software governed by the CeCILL license (see ./LICENSE)
61   !!----------------------------------------------------------------------
62CONTAINS
63
64#if defined key_iomput
65   !!----------------------------------------------------------------------
66   !!   'key_iomput'                                        use IOM library
67   !!----------------------------------------------------------------------
68   INTEGER FUNCTION dia_wri_alloc()
69      !
70      dia_wri_alloc = 0
71      !
72   END FUNCTION dia_wri_alloc
73
74
75   SUBROUTINE dia_wri( kt, Kmm )
76      !!---------------------------------------------------------------------
77      !!                  ***  ROUTINE dia_wri  ***
78      !!
79      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
80      !!      NETCDF format is used by default
81      !!     STATION_ASF
82      !!
83      !! ** Method  :  use iom_put
84      !!----------------------------------------------------------------------
85      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
86      INTEGER, INTENT( in ) ::   Kmm     ! ocean time level index
87      !!----------------------------------------------------------------------
88      !
89      IF( ln_timing )   CALL timing_start('dia_wri')
90      !
91      ! Output the initial state and forcings
92      IF( ninist == 1 ) THEN
93         CALL dia_wri_state( Kmm, 'output.init' )
94         ninist = 0
95      ENDIF
96      !
97      CALL iom_put(  "sst", sst_m(:,:) )    ! bulk surface temperature
98      CALL iom_put(  "sss", sss_m(:,:) )    ! surface salinity
99      !
100      CALL iom_put(  "ssu", ssu_m(:,:) )    ! ocean surface current along i-axis
101      CALL iom_put(  "ssv", ssv_m(:,:) )    ! ocean surface current along j-axis
102      !
103      IF( ln_timing )   CALL timing_stop('dia_wri')
104      !
105   END SUBROUTINE dia_wri
106
107#else
108   !!----------------------------------------------------------------------
109   !!   Default option                                  use IOIPSL  library
110   !!----------------------------------------------------------------------
111
112   INTEGER FUNCTION dia_wri_alloc()
113      !!----------------------------------------------------------------------
114      INTEGER, DIMENSION(2) :: ierr
115      !!----------------------------------------------------------------------
116      ierr = 0
117      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
118         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
119         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
120      !
121      dia_wri_alloc = MAXVAL(ierr)
122      CALL mpp_sum( 'diawri', dia_wri_alloc )
123      !
124   END FUNCTION dia_wri_alloc
125
126
127
128   SUBROUTINE dia_wri( kt, Kmm )
129      !!---------------------------------------------------------------------
130      !!                  ***  ROUTINE dia_wri  ***
131      !!
132      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
133      !!      NETCDF format is used by default
134      !!
135      !! ** Method  :   At the beginning of the first time step (nit000),
136      !!      define all the NETCDF files and fields
137      !!      At each time step call histdef to compute the mean if ncessary
138      !!      Each nn_write time step, output the instantaneous or mean fields
139      !!----------------------------------------------------------------------
140      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
141      INTEGER, INTENT( in ) ::   Kmm  ! ocean time level index
142      !
143      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
144      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
145      INTEGER  ::   inum = 11                                ! temporary logical unit
146      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
147      INTEGER  ::   ierr                                     ! error code return from allocation
148      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
149      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
150      !!----------------------------------------------------------------------
151      !
152      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==!
153         CALL dia_wri_state( Kmm, 'output.init' )
154         ninist = 0
155      ENDIF
156      !
157      IF( nn_write == -1 )   RETURN   ! we will never do any output
158      !
159      IF( ln_timing )   CALL timing_start('dia_wri')
160      !
161      ! 0. Initialisation
162      ! -----------------
163
164      ll_print = .FALSE.                  ! local variable for debugging
165      ll_print = ll_print .AND. lwp
166
167      ! Define frequency of output and means
168      clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes)
169#if defined key_diainstant
170      zsto = nn_write * rn_Dt
171      clop = "inst("//TRIM(clop)//")"
172#else
173      zsto=rn_Dt
174      clop = "ave("//TRIM(clop)//")"
175#endif
176      zout = nn_write * rn_Dt
177      zmax = ( nitend - nit000 + 1 ) * rn_Dt
178
179      ! Define indices of the horizontal output zoom and vertical limit storage
180      iimi = 1      ;      iima = jpi
181      ijmi = 1      ;      ijma = jpj
182      ipk = jpk
183
184      ! define time axis
185      it = kt
186      itmod = kt - nit000 + 1
187
188
189      ! 1. Define NETCDF files and fields at beginning of first time step
190      ! -----------------------------------------------------------------
191
192      IF( kt == nit000 ) THEN
193
194         ! Define the NETCDF files (one per grid)
195
196         ! Compute julian date from starting date of the run
197         CALL ymds2ju( nyear, nmonth, nday, rn_Dt, zjulian )
198         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
199         IF(lwp)WRITE(numout,*)
200         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
201            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
202         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
203            ' limit storage in depth = ', ipk
204
205         ! WRITE root name in date.file for use by postpro
206         IF(lwp) THEN
207            CALL dia_nam( clhstnam, nn_write,' ' )
208            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
209            WRITE(inum,*) clhstnam
210            CLOSE(inum)
211         ENDIF
212
213         ! Define the T grid FILE ( nid_T )
214
215         CALL dia_nam( clhstnam, nn_write, 'grid_T' )
216         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
217         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
218            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
219            &          nit000-1, zjulian, rn_Dt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
220         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
221            &           "m", ipk, gdept_1d, nz_T, "down" )
222         !                                                            ! Index of ocean points
223         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
224
225         ! Define the U grid FILE ( nid_U )
226
227         CALL dia_nam( clhstnam, nn_write, 'grid_U' )
228         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
229         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
230            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
231            &          nit000-1, zjulian, rn_Dt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
232         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
233            &           "m", ipk, gdept_1d, nz_U, "down" )
234         !                                                            ! Index of ocean points
235         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
236
237         ! Define the V grid FILE ( nid_V )
238
239         CALL dia_nam( clhstnam, nn_write, 'grid_V' )                   ! filename
240         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
241         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
242            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
243            &          nit000-1, zjulian, rn_Dt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
244         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
245            &          "m", ipk, gdept_1d, nz_V, "down" )
246         !                                                            ! Index of ocean points
247         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
248
249         ! No W grid FILE
250
251         ! Declare all the output fields as NETCDF variables
252
253         !                                                                                      !!! nid_T : 3D
254         CALL histdef( nid_T, "sst_m", "Sea Surface temperature"            , "C"      ,   &  ! sst
255            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
256         CALL histdef( nid_T, "sss_m", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
257            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
258         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
259            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
260         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! (sfx)
261            &         jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
262         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
263            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
264         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
265            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
266         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
267            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
268         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
269            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
270
271         CALL histend( nid_T, snc4chunks=snc4set )
272
273         !                                                                                      !!! nid_U : 3D
274         CALL histdef( nid_U, "ssu_m", "Velocity component in x-direction", "m/s"   ,         &  ! ssu
275            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
276         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
277            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
278
279         CALL histend( nid_U, snc4chunks=snc4set )
280
281         !                                                                                      !!! nid_V : 3D
282         CALL histdef( nid_V, "ssv_m", "Velocity component in y-direction", "m/s",            &  ! ssv_m
283            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
284         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
285            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
286
287         CALL histend( nid_V, snc4chunks=snc4set )
288
289         IF(lwp) WRITE(numout,*)
290         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
291         IF(ll_print) CALL FLUSH(numout )
292
293      ENDIF
294
295      ! 2. Start writing data
296      ! ---------------------
297
298      ! ndex(1) est utilise ssi l'avant dernier argument est different de
299      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
300      ! donne le nombre d'elements, et ndex la liste des indices a sortir
301
302      IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN
303         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
304         WRITE(numout,*) '~~~~~~ '
305      ENDIF
306
307      ! Write fields on T grid
308      CALL histwrite( nid_T, "sst_m", it, sst_m, ndim_hT, ndex_hT )   ! sea surface temperature
309      CALL histwrite( nid_T, "sss_m", it, sss_m, ndim_hT, ndex_hT )   ! sea surface salinity
310      CALL histwrite( nid_T, "sowaflup", it, (emp - rnf )  , ndim_hT, ndex_hT )   ! upward water flux
311      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
312      ! (includes virtual salt flux beneath ice
313      ! in linear free surface case)
314
315      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
316      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
317      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction
318      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed
319
320      ! Write fields on U grid
321      CALL histwrite( nid_U, "ssu_m"   , it, ssu_m         , ndim_hU, ndex_hU )   ! i-current speed
322      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
323
324      ! Write fields on V grid
325      CALL histwrite( nid_V, "ssv_m"   , it, ssv_m         , ndim_hV, ndex_hV )   ! j-current speed
326      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
327
328      ! 3. Close all files
329      ! ---------------------------------------
330      IF( kt == nitend ) THEN
331         CALL histclo( nid_T )
332         CALL histclo( nid_U )
333         CALL histclo( nid_V )
334      ENDIF
335      !
336      IF( ln_timing )   CALL timing_stop('dia_wri')
337      !
338   END SUBROUTINE dia_wri
339#endif
340
341   SUBROUTINE dia_wri_state( Kmm, cdfile_name )
342      !!---------------------------------------------------------------------
343      !!                 ***  ROUTINE dia_wri_state  ***
344      !!
345      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
346      !!      the instantaneous ocean state and forcing fields.
347      !!        Used to find errors in the initial state or save the last
348      !!      ocean state in case of abnormal end of a simulation
349      !!
350      !! ** Method  :   NetCDF files using ioipsl
351      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
352      !!      File 'output.abort.nc' is created in case of abnormal job end
353      !!----------------------------------------------------------------------
354      INTEGER           , INTENT( in ) ::   Kmm              ! time level index
355      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
356      !!
357      INTEGER :: inum
358      !!----------------------------------------------------------------------
359      !
360      IF(lwp) WRITE(numout,*)
361      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
362      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
363      IF(lwp) WRITE(numout,*) '                and named :', cdfile_name, '...nc'
364
365#if defined key_si3
366      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
367#else
368      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
369#endif
370
371      CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) )    ! now temperature
372      CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) )    ! now salinity
373      CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm)              )    ! sea surface height
374      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)                )    ! now i-velocity
375      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)                )    ! now j-velocity
376         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww             )    ! now k-velocity
377      CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf         )    ! freshwater budget
378      CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns         )    ! total heat flux
379      CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr               )    ! solar heat flux
380      CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i              )    ! ice fraction
381      CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
382      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
383
384#if defined key_si3
385      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid
386         CALL ice_wri_state( inum )
387      ENDIF
388#endif
389      !
390      CALL iom_close( inum )
391      !
392   END SUBROUTINE dia_wri_state
393
394   !!======================================================================
395END MODULE diawri
Note: See TracBrowser for help on using the repository browser.