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.
diawri.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/tests/STATION_ASF/MY_SRC – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/tests/STATION_ASF/MY_SRC/diawri.F90 @ 12254

Last change on this file since 12254 was 12249, checked in by laurent, 4 years ago

Made STATION_ASF testcase fully compliant with new timestepping scheme.

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