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 @ 11637

Last change on this file since 11637 was 11637, checked in by laurent, 5 years ago

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