source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/STATION_ASF/MY_SRC/diawri.F90 @ 13197

Last change on this file since 13197 was 13197, checked in by gsamson, 3 months ago

merge with trunk@r13136 with a more recent svn version; pass all SETTE tests; results identical to trunk@r13136; ticket #2419

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