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/src/SAS – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/SAS/diawri.F90 @ 11949

Last change on this file since 11949 was 11949, checked in by acc, 4 years ago

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