source: NEMO/releases/r4.0/r4.0-HEAD/src/SAS/diawri.F90 @ 12640

Last change on this file since 12640 was 11536, checked in by smasson, 21 months ago

trunk: merge dev_r10984_HPC-13 into the trunk

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