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

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/SAS/diawri.F90 @ 11463

Last change on this file since 11463 was 11053, checked in by davestorkey, 5 years ago

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