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

Last change on this file since 10425 was 10425, checked in by smasson, 2 years ago

trunk: merge back dev_r10164_HPC09_ESIWACE_PREP_MERGE@10424 into the trunk

  • Property svn:keywords set to Id
File size: 18.8 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 nwrite 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      IF( ln_timing )   CALL timing_start('dia_wri')
142      !
143      ! Output the initial state and forcings
144      IF( ninist == 1 ) THEN                       
145         CALL dia_wri_state( 'output.init' )
146         ninist = 0
147      ENDIF
148      !
149      ! 0. Initialisation
150      ! -----------------
151
152      ! local variable for debugging
153      ll_print = .FALSE.
154      ll_print = ll_print .AND. lwp
155
156      ! Define frequency of output and means
157      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
158      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cpu time)
159      ENDIF
160#if defined key_diainstant
161      zsto = nwrite * rdt
162      clop = "inst("//TRIM(clop)//")"
163#else
164      zsto=rdt
165      clop = "ave("//TRIM(clop)//")"
166#endif
167      zout = nwrite * rdt
168      zmax = ( nitend - nit000 + 1 ) * rdt
169
170      ! Define indices of the horizontal output zoom and vertical limit storage
171      iimi = 1      ;      iima = jpi
172      ijmi = 1      ;      ijma = jpj
173      ipk = jpk
174
175      ! define time axis
176      it = kt
177      itmod = kt - nit000 + 1
178
179
180      ! 1. Define NETCDF files and fields at beginning of first time step
181      ! -----------------------------------------------------------------
182
183      IF( kt == nit000 ) THEN
184
185         ! Define the NETCDF files (one per grid)
186
187         ! Compute julian date from starting date of the run
188         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
189         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
190         IF(lwp)WRITE(numout,*)
191         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
192            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
193         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
194                                 ' limit storage in depth = ', ipk
195
196         ! WRITE root name in date.file for use by postpro
197         IF(lwp) THEN
198            CALL dia_nam( clhstnam, nwrite,' ' )
199            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
200            WRITE(inum,*) clhstnam
201            CLOSE(inum)
202         ENDIF
203
204         ! Define the T grid FILE ( nid_T )
205
206         CALL dia_nam( clhstnam, nwrite, 'grid_T' )
207         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
208         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
209            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
210            &          nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
211         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
212            &           "m", ipk, gdept_1d, nz_T, "down" )
213         !                                                            ! Index of ocean points
214         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
215
216         ! Define the U grid FILE ( nid_U )
217
218         CALL dia_nam( clhstnam, nwrite, 'grid_U' )
219         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
220         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
221            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
222            &          nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
223         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
224            &           "m", ipk, gdept_1d, nz_U, "down" )
225         !                                                            ! Index of ocean points
226         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
227
228         ! Define the V grid FILE ( nid_V )
229
230         CALL dia_nam( clhstnam, nwrite, 'grid_V' )                   ! filename
231         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
232         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
233            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
234            &          nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
235         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
236            &          "m", ipk, gdept_1d, nz_V, "down" )
237         !                                                            ! Index of ocean points
238         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
239
240         ! No W grid FILE
241
242         ! Declare all the output fields as NETCDF variables
243
244         !                                                                                      !!! nid_T : 3D
245         CALL histdef( nid_T, "sst_m", "Sea Surface temperature"            , "C"      ,   &  ! sst
246            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
247         CALL histdef( nid_T, "sss_m", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
248            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
249         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
250            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
251         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! (sfx)
252             &         jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
253         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
254            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
255         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
256            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
257         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
258            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
259         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
260            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
261
262         CALL histend( nid_T, snc4chunks=snc4set )
263
264         !                                                                                      !!! nid_U : 3D
265         CALL histdef( nid_U, "ssu_m", "Velocity component in x-direction", "m/s"   ,         &  ! ssu
266            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
267         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
268            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
269
270         CALL histend( nid_U, snc4chunks=snc4set )
271
272         !                                                                                      !!! nid_V : 3D
273         CALL histdef( nid_V, "ssv_m", "Velocity component in y-direction", "m/s",            &  ! ssv_m
274            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
275         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
276            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
277
278         CALL histend( nid_V, snc4chunks=snc4set )
279
280         IF(lwp) WRITE(numout,*)
281         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
282         IF(ll_print) CALL FLUSH(numout )
283
284      ENDIF
285
286      ! 2. Start writing data
287      ! ---------------------
288
289      ! ndex(1) est utilise ssi l'avant dernier argument est diffferent de
290      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
291      ! donne le nombre d'elements, et ndex la liste des indices a sortir
292
293      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
294         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
295         WRITE(numout,*) '~~~~~~ '
296      ENDIF
297
298      ! Write fields on T grid
299      CALL histwrite( nid_T, "sst_m", it, sst_m, ndim_hT, ndex_hT )   ! sea surface temperature
300      CALL histwrite( nid_T, "sss_m", it, sss_m, ndim_hT, ndex_hT )   ! sea surface salinity
301      CALL histwrite( nid_T, "sowaflup", it, (emp - rnf )  , ndim_hT, ndex_hT )   ! upward water flux
302      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
303                                                                                  ! (includes virtual salt flux beneath ice
304                                                                                  ! in linear free surface case)
305
306      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
307      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
308      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
309      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
310
311         ! Write fields on U grid
312      CALL histwrite( nid_U, "ssu_m"   , it, ssu_m         , ndim_hU, ndex_hU )   ! i-current speed
313      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
314
315         ! Write fields on V grid
316      CALL histwrite( nid_V, "ssv_m"   , it, ssv_m         , ndim_hV, ndex_hV )   ! j-current speed
317      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
318
319      ! 3. Close all files
320      ! ---------------------------------------
321      IF( kt == nitend ) THEN
322         CALL histclo( nid_T )
323         CALL histclo( nid_U )
324         CALL histclo( nid_V )
325      ENDIF
326      !
327      IF( ln_timing )   CALL timing_stop('dia_wri')
328      !
329   END SUBROUTINE dia_wri
330#endif
331
332   SUBROUTINE dia_wri_state( cdfile_name )
333      !!---------------------------------------------------------------------
334      !!                 ***  ROUTINE dia_wri_state  ***
335      !!       
336      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
337      !!      the instantaneous ocean state and forcing fields.
338      !!        Used to find errors in the initial state or save the last
339      !!      ocean state in case of abnormal end of a simulation
340      !!
341      !! ** Method  :   NetCDF files using ioipsl
342      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
343      !!      File 'output.abort.nc' is created in case of abnormal job end
344      !!----------------------------------------------------------------------
345      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
346      !!
347      INTEGER :: inum
348      !!----------------------------------------------------------------------
349      !
350      IF(lwp) WRITE(numout,*)
351      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
352      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
353      IF(lwp) WRITE(numout,*) '                and named :', cdfile_name, '...nc'
354
355#if defined key_si3
356     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
357#else
358     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
359#endif
360
361      CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) )    ! now temperature
362      CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) )    ! now salinity
363      CALL iom_rstput( 0, 0, inum, 'sossheig', sshn              )    ! sea surface height
364      CALL iom_rstput( 0, 0, inum, 'vozocrtx', un                )    ! now i-velocity
365      CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity
366      CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn                )    ! now k-velocity
367      CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf         )    ! freshwater budget
368      CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns         )    ! total heat flux
369      CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr               )    ! solar heat flux
370      CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i              )    ! ice fraction
371      CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
372      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
373 
374#if defined key_si3
375      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid
376         CALL ice_wri_state( inum )
377      ENDIF
378#endif
379      !
380      CALL iom_close( inum )
381      !
382   END SUBROUTINE dia_wri_state
383
384   !!======================================================================
385END MODULE diawri
Note: See TracBrowser for help on using the repository browser.