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_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC – NEMO

source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/MY_SRC/diawri.F90 @ 11838

Last change on this file since 11838 was 11831, checked in by laurent, 4 years ago

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