source: NEMO/trunk/src/SAS/diawri.F90 @ 10068

Last change on this file since 10068 was 10068, checked in by nicolasmartin, 2 years ago

First part of modifications to have a common default header : fix typos and SVN keywords properties

  • Property svn:keywords set to Id
File size: 21.4 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   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   dia_wri       : create the standart output files
23   !!   dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
24   !!----------------------------------------------------------------------
25   USE oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE zdf_oce         ! ocean vertical physics
28   USE sbc_oce         ! Surface boundary condition: ocean fields
29   USE sbc_ice         ! Surface boundary condition: ice fields
30   USE sbcssr          ! restoring term toward SST/SSS climatology
31   USE phycst          ! physical constants
32   USE zdfmxl          ! mixed layer
33   USE dianam          ! build name of file (routine)
34   USE zdfddm          ! vertical  physics: double diffusion
35   USE diahth          ! thermocline diagnostics
36   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
37   USE in_out_manager  ! I/O manager
38   USE iom
39   USE ioipsl
40#if defined key_si3
41   USE icewri
42#endif
43   USE lib_mpp         ! MPP library
44   USE timing          ! preformance summary
45
46   IMPLICIT NONE
47   PRIVATE
48
49   PUBLIC   dia_wri                 ! routines called by step.F90
50   PUBLIC   dia_wri_state
51   PUBLIC   dia_wri_alloc           ! Called by nemogcm module
52
53   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
54   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
55   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
56   INTEGER ::   ndex(1)                              ! ???
57   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
58
59   !! * Substitutions
60#  include "vectopt_loop_substitute.h90"
61   !!----------------------------------------------------------------------
62   !! NEMO/SAS 4.0 , NEMO Consortium (2018)
63   !! $Id$
64   !! Software governed by the CeCILL license (see ./LICENSE)
65   !!----------------------------------------------------------------------
66CONTAINS
67
68# if defined key_iomput
69   !!----------------------------------------------------------------------
70   !!   'key_iomput'                                        use IOM library
71   !!----------------------------------------------------------------------
72   INTEGER FUNCTION dia_wri_alloc()
73      !
74      dia_wri_alloc = 0
75      !
76   END FUNCTION dia_wri_alloc
77
78   
79   SUBROUTINE dia_wri( kt )
80      !!---------------------------------------------------------------------
81      !!                  ***  ROUTINE dia_wri  ***
82      !!                   
83      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
84      !!      NETCDF format is used by default
85      !!      Standalone surface scheme
86      !!
87      !! ** Method  :  use iom_put
88      !!----------------------------------------------------------------------
89      !!
90      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
91      !!----------------------------------------------------------------------
92      !
93      ! Output the initial state and forcings
94      IF( ninist == 1 ) THEN
95         CALL dia_wri_state( 'output.init', kt )
96         ninist = 0
97      ENDIF
98      !
99   END SUBROUTINE dia_wri
100
101#else
102   !!----------------------------------------------------------------------
103   !!   Default option                                  use IOIPSL  library
104   !!----------------------------------------------------------------------
105   INTEGER FUNCTION dia_wri_alloc()
106      !!----------------------------------------------------------------------
107      INTEGER :: ierr
108      !!----------------------------------------------------------------------
109      !
110      ALLOCATE( ndex_hT(jpi*jpj), ndex_hU(jpi*jpj), ndex_hV(jpi*jpj), STAT=dia_wri_alloc )
111      IF( lk_mpp )   CALL mpp_sum( dia_wri_alloc )
112      !
113   END FUNCTION dia_wri_alloc
114   
115 
116   SUBROUTINE dia_wri( kt )
117      !!---------------------------------------------------------------------
118      !!                  ***  ROUTINE dia_wri  ***
119      !!                   
120      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
121      !!      NETCDF format is used by default
122      !!
123      !! ** Method  :   At the beginning of the first time step (nit000),
124      !!      define all the NETCDF files and fields
125      !!      At each time step call histdef to compute the mean if ncessary
126      !!      Each nwrite time step, output the instantaneous or mean fields
127      !!----------------------------------------------------------------------
128      !!
129      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
130      !!
131      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
132      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
133      INTEGER  ::   inum = 11                                ! temporary logical unit
134      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
135      INTEGER  ::   ierr                                     ! error code return from allocation
136      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
137      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
138      !!----------------------------------------------------------------------
139      !
140      IF( ln_timing )   CALL timing_start('dia_wri')
141      !
142      ! Output the initial state and forcings
143      IF( ninist == 1 ) THEN                       
144         CALL dia_wri_state( 'output.init', kt )
145         ninist = 0
146      ENDIF
147      !
148      ! 0. Initialisation
149      ! -----------------
150
151      ! local variable for debugging
152      ll_print = .FALSE.
153      ll_print = ll_print .AND. lwp
154
155      ! Define frequency of output and means
156      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
157      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cpu time)
158      ENDIF
159#if defined key_diainstant
160      zsto = nwrite * rdt
161      clop = "inst("//TRIM(clop)//")"
162#else
163      zsto=rdt
164      clop = "ave("//TRIM(clop)//")"
165#endif
166      zout = nwrite * rdt
167      zmax = ( nitend - nit000 + 1 ) * rdt
168
169      ! Define indices of the horizontal output zoom and vertical limit storage
170      iimi = 1      ;      iima = jpi
171      ijmi = 1      ;      ijma = jpj
172      ipk = jpk
173
174      ! define time axis
175      it = kt
176      itmod = kt - nit000 + 1
177
178
179      ! 1. Define NETCDF files and fields at beginning of first time step
180      ! -----------------------------------------------------------------
181
182      IF( kt == nit000 ) THEN
183
184         ! Define the NETCDF files (one per grid)
185
186         ! Compute julian date from starting date of the run
187         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
188         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
189         IF(lwp)WRITE(numout,*)
190         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
191            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
192         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
193                                 ' limit storage in depth = ', ipk
194
195         ! WRITE root name in date.file for use by postpro
196         IF(lwp) THEN
197            CALL dia_nam( clhstnam, nwrite,' ' )
198            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
199            WRITE(inum,*) clhstnam
200            CLOSE(inum)
201         ENDIF
202
203         ! Define the T grid FILE ( nid_T )
204
205         CALL dia_nam( clhstnam, nwrite, 'grid_T' )
206         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
207         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
208            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
209            &          nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
210         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
211            &           "m", ipk, gdept_1d, nz_T, "down" )
212         !                                                            ! Index of ocean points
213         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
214
215         ! Define the U grid FILE ( nid_U )
216
217         CALL dia_nam( clhstnam, nwrite, 'grid_U' )
218         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
219         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
220            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
221            &          nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
222         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
223            &           "m", ipk, gdept_1d, nz_U, "down" )
224         !                                                            ! Index of ocean points
225         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
226
227         ! Define the V grid FILE ( nid_V )
228
229         CALL dia_nam( clhstnam, nwrite, 'grid_V' )                   ! filename
230         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
231         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
232            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
233            &          nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
234         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
235            &          "m", ipk, gdept_1d, nz_V, "down" )
236         !                                                            ! Index of ocean points
237         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
238
239         ! No W grid FILE
240
241         ! Declare all the output fields as NETCDF variables
242
243         !                                                                                      !!! nid_T : 3D
244         CALL histdef( nid_T, "sst_m", "Sea Surface temperature"            , "C"      ,   &  ! sst
245            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
246         CALL histdef( nid_T, "sss_m", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
247            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
248         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
249            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
250         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! (sfx)
251             &         jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
252         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
253            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
254         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
255            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
256         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
257            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
258         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
259            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
260
261         CALL histend( nid_T, snc4chunks=snc4set )
262
263         !                                                                                      !!! nid_U : 3D
264         CALL histdef( nid_U, "ssu_m", "Velocity component in x-direction", "m/s"   ,         &  ! ssu
265            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
266         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
267            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
268
269         CALL histend( nid_U, snc4chunks=snc4set )
270
271         !                                                                                      !!! nid_V : 3D
272         CALL histdef( nid_V, "ssv_m", "Velocity component in y-direction", "m/s",            &  ! ssv_m
273            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
274         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
275            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
276
277         CALL histend( nid_V, snc4chunks=snc4set )
278
279         IF(lwp) WRITE(numout,*)
280         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
281         IF(ll_print) CALL FLUSH(numout )
282
283      ENDIF
284
285      ! 2. Start writing data
286      ! ---------------------
287
288      ! ndex(1) est utilise ssi l'avant dernier argument est diffferent de
289      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
290      ! donne le nombre d'elements, et ndex la liste des indices a sortir
291
292      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
293         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
294         WRITE(numout,*) '~~~~~~ '
295      ENDIF
296
297      ! Write fields on T grid
298      CALL histwrite( nid_T, "sst_m", it, sst_m, ndim_hT, ndex_hT )   ! sea surface temperature
299      CALL histwrite( nid_T, "sss_m", it, sss_m, ndim_hT, ndex_hT )   ! sea surface salinity
300      CALL histwrite( nid_T, "sowaflup", it, (emp - rnf )  , ndim_hT, ndex_hT )   ! upward water flux
301      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
302                                                                                  ! (includes virtual salt flux beneath ice
303                                                                                  ! in linear free surface case)
304
305      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
306      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
307      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
308      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
309
310         ! Write fields on U grid
311      CALL histwrite( nid_U, "ssu_m"   , it, ssu_m         , ndim_hU, ndex_hU )   ! i-current speed
312      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
313
314         ! Write fields on V grid
315      CALL histwrite( nid_V, "ssv_m"   , it, ssv_m         , ndim_hV, ndex_hV )   ! j-current speed
316      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
317
318      ! 3. Close all files
319      ! ---------------------------------------
320      IF( kt == nitend ) THEN
321         CALL histclo( nid_T )
322         CALL histclo( nid_U )
323         CALL histclo( nid_V )
324      ENDIF
325      !
326      IF( ln_timing )   CALL timing_stop('dia_wri')
327      !
328   END SUBROUTINE dia_wri
329#endif
330
331   SUBROUTINE dia_wri_state( cdfile_name, kt )
332      !!---------------------------------------------------------------------
333      !!                 ***  ROUTINE dia_wri_state  ***
334      !!       
335      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
336      !!      the instantaneous ocean state and forcing fields.
337      !!        Used to find errors in the initial state or save the last
338      !!      ocean state in case of abnormal end of a simulation
339      !!
340      !! ** Method  :   NetCDF files using ioipsl
341      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
342      !!      File 'output.abort.nc' is created in case of abnormal job end
343      !!----------------------------------------------------------------------
344      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
345      INTEGER           , INTENT( in ) ::   kt               ! ocean time-step index
346      !!
347      CHARACTER (len=32) :: clname
348      CHARACTER (len=40) :: clop
349      INTEGER  ::   id_i , nz_i, nh_i       
350      INTEGER, DIMENSION(1) ::   idex             ! local workspace
351      REAL(wp) ::   zsto, zout, zmax, zjulian
352      !!----------------------------------------------------------------------
353      !
354      IF( ln_timing )   CALL timing_start('dia_wri_state')
355
356      ! 0. Initialisation
357      ! -----------------
358
359      ! Define name, frequency of output and means
360      clname = cdfile_name
361      IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
362      zsto = rdt
363      clop = "inst(x)"           ! no use of the mask value (require less cpu time)
364      zout = rdt
365      zmax = ( nitend - nit000 + 1 ) * rdt
366
367      IF(lwp) WRITE(numout,*)
368      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
369      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
370      IF(lwp) WRITE(numout,*) '                and named :', clname, '.nc'
371
372
373      ! 1. Define NETCDF files and fields at beginning of first time step
374      ! -----------------------------------------------------------------
375
376      ! Compute julian date from starting date of the run
377      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )         ! time axis
378      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
379      CALL histbeg( clname, jpi, glamt, jpj, gphit,   &
380          1, jpi, 1, jpj, nit000-1, zjulian, rdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
381      CALL histvert( id_i, "deptht", "Vertical T levels",   &    ! Vertical grid : gdept
382          "m", jpk, gdept_1d, nz_i, "down")
383
384      ! Declare all the output fields as NetCDF variables
385
386      CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S",   &   ! net freshwater
387         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
388      CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2"   ,   &   ! net heat flux
389         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
390      CALL histdef( id_i, "soshfldo", "Shortwave Radiation"   , "W/m2"   ,   &   ! solar flux
391         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
392      CALL histdef( id_i, "soicecov", "Ice fraction"          , "[0,1]"  ,   &   ! fr_i
393         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
394      CALL histdef( id_i, "sozotaux", "Zonal Wind Stress"     , "N/m2"   ,   &   ! i-wind stress
395         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
396      CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2"   ,   &   ! j-wind stress
397         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
398
399#if defined key_si3
400      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + lim but no-ice in child grid
401         CALL ice_wri_state( kt, id_i, nh_i )
402      ENDIF
403#else
404      CALL histend( id_i, snc4chunks=snc4set )
405#endif
406
407      ! 2. Start writing data
408      ! ---------------------
409      ! idex(1) est utilise ssi l'avant dernier argument est diffferent de
410      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
411      ! donne le nombre d'elements, et idex la liste des indices a sortir
412      idex(1) = 1   ! init to avoid compil warning
413
414      ! Write all fields on T grid
415      CALL histwrite( id_i, "sowaflup", kt, emp              , jpi*jpj    , idex )    ! freshwater budget
416      CALL histwrite( id_i, "sohefldo", kt, qsr + qns        , jpi*jpj    , idex )    ! total heat flux
417      CALL histwrite( id_i, "soshfldo", kt, qsr              , jpi*jpj    , idex )    ! solar heat flux
418      CALL histwrite( id_i, "soicecov", kt, fr_i             , jpi*jpj    , idex )    ! ice fraction
419      CALL histwrite( id_i, "sozotaux", kt, utau             , jpi*jpj    , idex )    ! i-wind stress
420      CALL histwrite( id_i, "sometauy", kt, vtau             , jpi*jpj    , idex )    ! j-wind stress
421
422      ! 3. Close the file
423      ! -----------------
424      CALL histclo( id_i )
425#if ! defined key_iomput
426      IF( ninist /= 1  ) THEN
427         CALL histclo( nid_T )
428         CALL histclo( nid_U )
429         CALL histclo( nid_V )
430      ENDIF
431#endif
432      !
433      IF( ln_timing )   CALL timing_stop('dia_wri_state')
434      !
435   END SUBROUTINE dia_wri_state
436
437   !!======================================================================
438END MODULE diawri
Note: See TracBrowser for help on using the repository browser.