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 branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 2620

Last change on this file since 2620 was 2613, checked in by gm, 13 years ago

dynamic mem: #785 ; move the allocation of ice in iceini_2/iceini module + bug fixes (define key_esopa)
  • Property svn:keywords set to Id
File size: 39.3 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 ldftra_oce      ! ocean active tracers: lateral physics
29   USE ldfdyn_oce      ! ocean dynamics: lateral physics
30   USE sol_oce         ! solver variables
31   USE sbc_oce         ! Surface boundary condition: ocean fields
32   USE sbc_ice         ! Surface boundary condition: ice fields
33   USE sbcssr          ! restoring term toward SST/SSS climatology
34   USE phycst          ! physical constants
35   USE zdfmxl          ! mixed layer
36   USE dianam          ! build name of file (routine)
37   USE zdfddm          ! vertical  physics: double diffusion
38   USE diahth          ! thermocline diagnostics
39   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
40   USE in_out_manager  ! I/O manager
41   USE diadimg         ! dimg direct access file format output
42   USE diaar5, ONLY :   lk_diaar5
43   USE iom
44   USE ioipsl
45#if defined key_lim2
46   USE limwri_2 
47#endif
48   USE dtatem
49   USE dtasal
50
51   IMPLICIT NONE
52   PRIVATE
53
54   PUBLIC   dia_wri                 ! routines called by step.F90
55   PUBLIC   dia_wri_state
56   PUBLIC   dia_wri_alloc           ! Called by nemogcm module
57
58   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
59   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
60   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
61   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
62   INTEGER ::   ndex(1)                              ! ???
63   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
64   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
65
66   !! * Substitutions
67#  include "zdfddm_substitute.h90"
68#  include "domzgr_substitute.h90"
69#  include "vectopt_loop_substitute.h90"
70   !!----------------------------------------------------------------------
71   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
72   !! $Id $
73   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
74   !!----------------------------------------------------------------------
75CONTAINS
76
77   FUNCTION dia_wri_alloc()
78      !!----------------------------------------------------------------------
79      IMPLICIT none
80      INTEGER :: dia_wri_alloc
81      INTEGER, DIMENSION(2) :: ierr
82      !!----------------------------------------------------------------------
83      !
84      ierr = 0
85      !
86      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
87         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
88         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
89         !
90      dia_wri_alloc = MAXVAL(ierr)
91      IF( lk_mpp )   CALL mpp_sum( ierr )
92      !
93  END FUNCTION dia_wri_alloc
94
95#if defined key_dimgout
96   !!----------------------------------------------------------------------
97   !!   'key_dimgout'                                      DIMG output file
98   !!----------------------------------------------------------------------
99#   include "diawri_dimg.h90"
100
101#else
102   !!----------------------------------------------------------------------
103   !!   Default option                                   NetCDF output file
104   !!----------------------------------------------------------------------
105# if defined key_iomput
106   !!----------------------------------------------------------------------
107   !!   'key_iomput'                                        use IOM library
108   !!----------------------------------------------------------------------
109
110   SUBROUTINE dia_wri( kt )
111      !!---------------------------------------------------------------------
112      !!                  ***  ROUTINE dia_wri  ***
113      !!                   
114      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
115      !!      NETCDF format is used by default
116      !!
117      !! ** Method  :  use iom_put
118      !!----------------------------------------------------------------------
119      USE oce, ONLY :   z3d => ta   ! use ta as 3D workspace
120      USE wrk_nemo, ONLY: wrk_use, wrk_release
121      USE wrk_nemo, ONLY: z2d => wrk_2d_1
122      !!
123      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
124      !!
125      INTEGER                      ::   ji, jj, jk              ! dummy loop indices
126      REAL(wp)                     ::   zztmp, zztmpx, zztmpy   !
127      !!----------------------------------------------------------------------
128      !
129      IF( .not. wrk_use(2, 1))THEN
130         CALL ctl_stop('dia_wri: ERROR - requested 2D workspace unavailable.')
131         RETURN
132      END IF
133      !
134      ! Output the initial state and forcings
135      IF( ninist == 1 ) THEN                       
136         CALL dia_wri_state( 'output.init', kt )
137         ninist = 0
138      ENDIF
139
140      CALL iom_put( "toce"   , tn                    )    ! temperature
141      CALL iom_put( "soce"   , sn                    )    ! salinity
142      CALL iom_put( "sst"    , tn(:,:,1)             )    ! sea surface temperature
143      CALL iom_put( "sst2"   , tn(:,:,1) * tn(:,:,1) )    ! square of sea surface temperature
144      CALL iom_put( "sss"    , sn(:,:,1)             )    ! sea surface salinity
145      CALL iom_put( "sss2"   , sn(:,:,1) * sn(:,:,1) )    ! square of sea surface salinity
146      CALL iom_put( "uoce"   , un                    )    ! i-current     
147      CALL iom_put( "voce"   , vn                    )    ! j-current
148     
149      CALL iom_put( "avt"    , avt                   )    ! T vert. eddy diff. coef.
150      CALL iom_put( "avm"    , avmu                  )    ! T vert. eddy visc. coef.
151      IF( lk_zdfddm ) THEN
152         CALL iom_put( "avs" , fsavs(:,:,:)          )    ! S vert. eddy diff. coef.
153      ENDIF
154
155      DO jj = 2, jpjm1                                    ! sst gradient
156         DO ji = fs_2, fs_jpim1   ! vector opt.
157            zztmp      = tn(ji,jj,1)
158            zztmpx     = ( tn(ji+1,jj  ,1) - zztmp ) / e1u(ji,jj) + ( zztmp - tn(ji-1,jj  ,1) ) / e1u(ji-1,jj  )
159            zztmpy     = ( tn(ji  ,jj+1,1) - zztmp ) / e2v(ji,jj) + ( zztmp - tn(ji  ,jj-1,1) ) / e2v(ji  ,jj-1)
160            z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
161               &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
162         END DO
163      END DO
164      CALL lbc_lnk( z2d, 'T', 1. )
165      CALL iom_put( "sstgrad2",  z2d               )    ! square of module of sst gradient
166!CDIR NOVERRCHK
167      z2d(:,:) = SQRT( z2d(:,:) )
168      CALL iom_put( "sstgrad" ,  z2d               )    ! module of sst gradient
169
170      IF( lk_diaar5 ) THEN
171         z3d(:,:,jpk) = 0.e0
172         DO jk = 1, jpkm1
173            z3d(:,:,jk) = rau0 * un(:,:,jk) * e1u(:,:) * fse3u(:,:,jk)
174         END DO
175         CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction
176         zztmp = 0.5 * rcp
177         z2d(:,:) = 0.e0 
178         DO jk = 1, jpkm1
179            DO jj = 2, jpjm1
180               DO ji = fs_2, fs_jpim1   ! vector opt.
181                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * zztmp * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) )
182               END DO
183            END DO
184         END DO
185         CALL lbc_lnk( z2d, 'U', -1. )
186         CALL iom_put( "u_heattr", z2d )                  ! heat transport in i-direction
187         DO jk = 1, jpkm1
188            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e2v(:,:) * fse3v(:,:,jk)
189         END DO
190         CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction
191         z2d(:,:) = 0.e0 
192         DO jk = 1, jpkm1
193            DO jj = 2, jpjm1
194               DO ji = fs_2, fs_jpim1   ! vector opt.
195                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * zztmp * ( tn(ji,jj,jk) + tn(ji,jj+1,jk) )
196               END DO
197            END DO
198         END DO
199         CALL lbc_lnk( z2d, 'V', -1. )
200         CALL iom_put( "v_heattr", z2d )                  !  heat transport in i-direction
201      ENDIF
202      !
203      IF( .not. wrk_release(2, 1))THEN
204         CALL ctl_stop('dia_wri: ERROR - failed to release 2D workspace.')
205         RETURN
206      END IF
207      !
208   END SUBROUTINE dia_wri
209
210#else
211   !!----------------------------------------------------------------------
212   !!   Default option                                  use IOIPSL  library
213   !!----------------------------------------------------------------------
214
215   SUBROUTINE dia_wri( kt )
216      !!---------------------------------------------------------------------
217      !!                  ***  ROUTINE dia_wri  ***
218      !!                   
219      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
220      !!      NETCDF format is used by default
221      !!
222      !! ** Method  :   At the beginning of the first time step (nit000),
223      !!      define all the NETCDF files and fields
224      !!      At each time step call histdef to compute the mean if ncessary
225      !!      Each nwrite time step, output the instantaneous or mean fields
226      !!----------------------------------------------------------------------
227      USE wrk_nemo, ONLY: wrk_use, wrk_release
228      USE wrk_nemo, ONLY: zw2d => wrk_2d_1
229      !!
230      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
231      !!
232      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
233      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
234      INTEGER  ::   inum = 11                                ! temporary logical unit
235      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
236      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt           ! local scalars
237      !!----------------------------------------------------------------------
238      !
239      IF( .not. wrk_use(2, 1))THEN
240         CALL ctl_stop('dia_wri: ERROR - requested 2D workspace unavailable.')
241         RETURN
242      END IF
243      !
244      ! Output the initial state and forcings
245      IF( ninist == 1 ) THEN                       
246         CALL dia_wri_state( 'output.init', kt )
247         ninist = 0
248      ENDIF
249      !
250      ! 0. Initialisation
251      ! -----------------
252
253      ! local variable for debugging
254      ll_print = .FALSE.
255      ll_print = ll_print .AND. lwp
256
257      ! Define frequency of output and means
258      zdt = rdt
259      IF( nacc == 1 ) zdt = rdtmin
260      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
261      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cpu time)
262      ENDIF
263#if defined key_diainstant
264      zsto = nwrite * zdt
265      clop = "inst("//TRIM(clop)//")"
266#else
267      zsto=zdt
268      clop = "ave("//TRIM(clop)//")"
269#endif
270      zout = nwrite * zdt
271      zmax = ( nitend - nit000 + 1 ) * zdt
272
273      ! Define indices of the horizontal output zoom and vertical limit storage
274      iimi = 1      ;      iima = jpi
275      ijmi = 1      ;      ijma = jpj
276      ipk = jpk
277
278      ! define time axis
279      it = kt
280      itmod = kt - nit000 + 1
281
282
283      ! 1. Define NETCDF files and fields at beginning of first time step
284      ! -----------------------------------------------------------------
285
286      IF( kt == nit000 ) THEN
287
288         ! Define the NETCDF files (one per grid)
289
290         ! Compute julian date from starting date of the run
291         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
292         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
293         IF(lwp)WRITE(numout,*)
294         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
295            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
296         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
297                                 ' limit storage in depth = ', ipk
298
299         ! WRITE root name in date.file for use by postpro
300         IF(lwp) THEN
301            CALL dia_nam( clhstnam, nwrite,' ' )
302            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
303            WRITE(inum,*) clhstnam
304            CLOSE(inum)
305         ENDIF
306
307         ! Define the T grid FILE ( nid_T )
308
309         CALL dia_nam( clhstnam, nwrite, 'grid_T' )
310         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
311         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
312            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
313            &          nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
314         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
315            &           "m", ipk, gdept_0, nz_T, "down" )
316         !                                                            ! Index of ocean points
317         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
318         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
319
320         ! Define the U grid FILE ( nid_U )
321
322         CALL dia_nam( clhstnam, nwrite, 'grid_U' )
323         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
324         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
325            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
326            &          nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
327         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
328            &           "m", ipk, gdept_0, nz_U, "down" )
329         !                                                            ! Index of ocean points
330         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
331         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
332
333         ! Define the V grid FILE ( nid_V )
334
335         CALL dia_nam( clhstnam, nwrite, 'grid_V' )                   ! filename
336         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
337         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
338            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
339            &          nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
340         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
341            &          "m", ipk, gdept_0, nz_V, "down" )
342         !                                                            ! Index of ocean points
343         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
344         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
345
346         ! Define the W grid FILE ( nid_W )
347
348         CALL dia_nam( clhstnam, nwrite, 'grid_W' )                   ! filename
349         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
350         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
351            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
352            &          nit000-1, zjulian, zdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
353         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
354            &          "m", ipk, gdepw_0, nz_W, "down" )
355
356
357         ! Declare all the output fields as NETCDF variables
358
359         !                                                                                      !!! nid_T : 3D
360         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
361            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
362         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
363            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
364         !                                                                                      !!! nid_T : 2D
365         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
366            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
367         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
368            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
369         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
370            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
371!!$#if defined key_lim3 || defined key_lim2
372!!$         ! sowaflup = sowaflep + sorunoff + sowafldp + a term associated to
373!!$         !    internal damping to Levitus that can be diagnosed from others
374!!$         ! sowaflcd = sowaflep + sorunoff + sowafldp + iowaflup
375!!$         CALL histdef( nid_T, "iowaflup", "Ice=>ocean net freshwater"          , "kg/m2/s",   &  ! fsalt
376!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
377!!$         CALL histdef( nid_T, "sowaflep", "atmos=>ocean net freshwater"        , "kg/m2/s",   &  ! fmass
378!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
379!!$#endif
380         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
381            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
382!!$         CALL histdef( nid_T, "sorunoff", "Runoffs"                            , "Kg/m2/s",   &  ! runoffs
383!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
384         CALL histdef( nid_T, "sowaflcd", "concentration/dilution water flux"  , "kg/m2/s",   &  ! (emps-rnf)
385            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
386         CALL histdef( nid_T, "sosalflx", "Surface Salt Flux"                  , "Kg/m2/s",   &  ! (emps-rnf) * sn
387            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
388         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
389            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
390         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
391            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
392         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
393            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
394         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
395            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
396         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
397            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
398         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
399            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
400#if ! defined key_coupled
401         CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
402            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
403         CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
404            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
405         CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
406            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
407#endif
408
409
410
411#if ( defined key_coupled && ! defined key_lim3 && ! defined key_lim2 )
412         CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
413            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
414         CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
415            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
416         CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn
417            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
418#endif
419         clmx ="l_max(only(x))"    ! max index on a period
420         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
421            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
422#if defined key_diahth
423         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
424            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
425         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
426            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
427         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
428            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
429         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "W"      ,   & ! htc3
430            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
431#endif
432
433#if defined key_coupled 
434# if defined key_lim3
435         Must be adapted to LIM3
436# else
437         CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice
438            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
439         CALL histdef( nid_T,"soicealb" , "Ice Albedo"                         , "[0,1]"  ,   &  ! alb_ice
440            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
441# endif 
442#endif
443
444         CALL histend( nid_T, snc4chunks=snc4set )
445
446         !                                                                                      !!! nid_U : 3D
447         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
448            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
449#if defined key_diaeiv
450         CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
451            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
452#endif
453         !                                                                                      !!! nid_U : 2D
454         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
455            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
456
457         CALL histend( nid_U, snc4chunks=snc4set )
458
459         !                                                                                      !!! nid_V : 3D
460         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
461            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
462#if defined key_diaeiv
463         CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
464            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
465#endif
466         !                                                                                      !!! nid_V : 2D
467         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
468            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
469
470         CALL histend( nid_V, snc4chunks=snc4set )
471
472         !                                                                                      !!! nid_W : 3D
473         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
474            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
475#if defined key_diaeiv
476         CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
477            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
478#endif
479         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
480            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
481         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avmu
482            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
483
484         IF( lk_zdfddm ) THEN
485            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
486               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
487         ENDIF
488         !                                                                                      !!! nid_W : 2D
489#if defined key_traldf_c2d
490         CALL histdef( nid_W, "soleahtw", "lateral eddy diffusivity"           , "m2/s"   ,   &  ! ahtw
491            &          jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
492# if defined key_traldf_eiv 
493            CALL histdef( nid_W, "soleaeiw", "eddy induced vel. coeff. at w-point", "m2/s",   &  ! aeiw
494               &       jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
495# endif
496#endif
497
498         CALL histend( nid_W, snc4chunks=snc4set )
499
500         IF(lwp) WRITE(numout,*)
501         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
502         IF(ll_print) CALL FLUSH(numout )
503
504      ENDIF
505
506      ! 2. Start writing data
507      ! ---------------------
508
509      ! ndex(1) est utilise ssi l'avant dernier argument est diffferent de
510      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
511      ! donne le nombre d'elements, et ndex la liste des indices a sortir
512
513      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
514         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
515         WRITE(numout,*) '~~~~~~ '
516      ENDIF
517
518      ! Write fields on T grid
519      CALL histwrite( nid_T, "votemper", it, tn            , ndim_T , ndex_T  )   ! temperature
520      CALL histwrite( nid_T, "vosaline", it, sn            , ndim_T , ndex_T  )   ! salinity
521      CALL histwrite( nid_T, "sosstsst", it, tn(:,:,1)     , ndim_hT, ndex_hT )   ! sea surface temperature
522      CALL histwrite( nid_T, "sosaline", it, sn(:,:,1)     , ndim_hT, ndex_hT )   ! sea surface salinity
523      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
524!!$#if  defined key_lim3 || defined key_lim2
525!!$      CALL histwrite( nid_T, "iowaflup", it, fsalt(:,:)    , ndim_hT, ndex_hT )   ! ice=>ocean water flux
526!!$      CALL histwrite( nid_T, "sowaflep", it, fmass(:,:)    , ndim_hT, ndex_hT )   ! atmos=>ocean water flux
527!!$#endif
528      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
529!!$      CALL histwrite( nid_T, "sorunoff", it, runoff        , ndim_hT, ndex_hT )   ! runoff
530      CALL histwrite( nid_T, "sowaflcd", it, ( emps-rnf )  , ndim_hT, ndex_hT )   ! c/d water flux
531      zw2d(:,:) = ( emps(:,:) - rnf(:,:) ) * sn(:,:,1) * tmask(:,:,1)
532      CALL histwrite( nid_T, "sosalflx", it, zw2d          , ndim_hT, ndex_hT )   ! c/d salt flux
533      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
534      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
535      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
536      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
537      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
538      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
539#if ! defined key_coupled
540      CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
541      CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
542      zw2d(:,:) = erp(:,:) * sn(:,:,1) * tmask(:,:,1)
543      CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
544#endif
545#if ( defined key_coupled && ! defined key_lim3 && ! defined key_lim2 )
546      CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
547      CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
548         zw2d(:,:) = erp(:,:) * sn(:,:,1) * tmask(:,:,1)
549      CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
550#endif
551      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
552      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
553
554#if defined key_diahth
555      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
556      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
557      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
558      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
559#endif
560
561#if defined key_coupled 
562# if defined key_lim3
563      Must be adapted for LIM3
564      CALL histwrite( nid_T, "soicetem", it, tn_ice        , ndim_hT, ndex_hT )   ! surf. ice temperature
565      CALL histwrite( nid_T, "soicealb", it, alb_ice       , ndim_hT, ndex_hT )   ! ice albedo
566# else
567      CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature
568      CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo
569# endif
570#endif
571         ! Write fields on U grid
572      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
573#if defined key_diaeiv
574      CALL histwrite( nid_U, "vozoeivu", it, u_eiv         , ndim_U , ndex_U )    ! i-eiv current
575#endif
576      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
577
578         ! Write fields on V grid
579      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
580#if defined key_diaeiv
581      CALL histwrite( nid_V, "vomeeivv", it, v_eiv         , ndim_V , ndex_V  )   ! j-eiv current
582#endif
583      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
584
585         ! Write fields on W grid
586      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
587#   if defined key_diaeiv
588      CALL histwrite( nid_W, "voveeivw", it, w_eiv          , ndim_T, ndex_T )    ! vert. eiv current
589#   endif
590      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
591      CALL histwrite( nid_W, "votkeavm", it, avmu           , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
592      IF( lk_zdfddm ) THEN
593         CALL histwrite( nid_W, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T )    ! S vert. eddy diff. coef.
594      ENDIF
595#if defined key_traldf_c2d
596      CALL histwrite( nid_W, "soleahtw", it, ahtw          , ndim_hT, ndex_hT )   ! lateral eddy diff. coef.
597# if defined key_traldf_eiv
598         CALL histwrite( nid_W, "soleaeiw", it, aeiw       , ndim_hT, ndex_hT )   ! EIV coefficient at w-point
599# endif
600#endif
601
602      ! 3. Close all files
603      ! ---------------------------------------
604      IF( kt == nitend ) THEN
605         CALL histclo( nid_T )
606         CALL histclo( nid_U )
607         CALL histclo( nid_V )
608         CALL histclo( nid_W )
609      ENDIF
610      !
611      IF( .not. wrk_release(2, 1))THEN
612         CALL ctl_stop('dia_wri: ERROR - failed to release 2D workspace.')
613         RETURN
614      END IF
615      !
616   END SUBROUTINE dia_wri
617# endif
618
619#endif
620
621   SUBROUTINE dia_wri_state( cdfile_name, kt )
622      !!---------------------------------------------------------------------
623      !!                 ***  ROUTINE dia_wri_state  ***
624      !!       
625      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
626      !!      the instantaneous ocean state and forcing fields.
627      !!        Used to find errors in the initial state or save the last
628      !!      ocean state in case of abnormal end of a simulation
629      !!
630      !! ** Method  :   NetCDF files using ioipsl
631      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
632      !!      File 'output.abort.nc' is created in case of abnormal job end
633      !!----------------------------------------------------------------------
634      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
635      INTEGER           , INTENT( in ) ::   kt               ! ocean time-step index
636      !!
637      CHARACTER (len=32) :: clname
638      CHARACTER (len=40) :: clop
639      INTEGER  ::   id_i , nz_i, nh_i       
640      INTEGER, DIMENSION(1) ::   idex             ! local workspace
641      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt
642      !!----------------------------------------------------------------------
643
644      ! 0. Initialisation
645      ! -----------------
646
647      ! Define name, frequency of output and means
648      clname = cdfile_name
649      IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
650      zdt  = rdt
651      zsto = rdt
652      clop = "inst(x)"           ! no use of the mask value (require less cpu time)
653      zout = rdt
654      zmax = ( nitend - nit000 + 1 ) * zdt
655
656      IF(lwp) WRITE(numout,*)
657      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
658      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
659      IF(lwp) WRITE(numout,*) '                and named :', clname, '.nc'
660
661
662      ! 1. Define NETCDF files and fields at beginning of first time step
663      ! -----------------------------------------------------------------
664
665      ! Compute julian date from starting date of the run
666      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )         ! time axis
667      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
668      CALL histbeg( clname, jpi, glamt, jpj, gphit,   &
669          1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
670      CALL histvert( id_i, "deptht", "Vertical T levels",   &    ! Vertical grid : gdept
671          "m", jpk, gdept_0, nz_i, "down")
672
673      ! Declare all the output fields as NetCDF variables
674
675      CALL histdef( id_i, "vosaline", "Salinity"              , "PSU"    ,   &   ! salinity
676         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
677      CALL histdef( id_i, "votemper", "Temperature"           , "C"      ,   &   ! temperature
678         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
679      CALL histdef( id_i, "sossheig", "Sea Surface Height"    , "m"      ,   &  ! ssh
680         &          jpi, jpj, nh_i, 1  , 1, 1  , nz_i, 32, clop, zsto, zout )
681      CALL histdef( id_i, "vozocrtx", "Zonal Current"         , "m/s"    ,   &   ! zonal current
682         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
683      CALL histdef( id_i, "vomecrty", "Meridional Current"    , "m/s"    ,   &   ! meridonal current
684         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
685      CALL histdef( id_i, "vovecrtz", "Vertical Velocity"     , "m/s"    ,   &   ! vertical current
686         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
687      CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S",   &   ! net freshwater
688         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
689      CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2"   ,   &   ! net heat flux
690         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
691      CALL histdef( id_i, "soshfldo", "Shortwave Radiation"   , "W/m2"   ,   &   ! solar flux
692         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
693      CALL histdef( id_i, "soicecov", "Ice fraction"          , "[0,1]"  ,   &   ! fr_i
694         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
695      CALL histdef( id_i, "sozotaux", "Zonal Wind Stress"     , "N/m2"   ,   &   ! i-wind stress
696         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
697      CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2"   ,   &   ! j-wind stress
698         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
699
700#if defined key_lim2
701      CALL lim_wri_state_2( kt, id_i, nh_i )
702#else
703      CALL histend( id_i, snc4chunks=snc4set )
704#endif
705
706      ! 2. Start writing data
707      ! ---------------------
708      ! idex(1) est utilise ssi l'avant dernier argument est diffferent de
709      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
710      ! donne le nombre d'elements, et idex la liste des indices a sortir
711      idex(1) = 1   ! init to avoid compil warning
712
713      ! Write all fields on T grid
714      CALL histwrite( id_i, "votemper", kt, tn      , jpi*jpj*jpk, idex )    ! now temperature
715      CALL histwrite( id_i, "vosaline", kt, sn      , jpi*jpj*jpk, idex )    ! now salinity
716      CALL histwrite( id_i, "sossheig", kt, sshn     , jpi*jpj    , idex )    ! sea surface height
717      CALL histwrite( id_i, "vozocrtx", kt, un       , jpi*jpj*jpk, idex )    ! now i-velocity
718      CALL histwrite( id_i, "vomecrty", kt, vn       , jpi*jpj*jpk, idex )    ! now j-velocity
719      CALL histwrite( id_i, "vovecrtz", kt, wn       , jpi*jpj*jpk, idex )    ! now k-velocity
720      CALL histwrite( id_i, "sowaflup", kt, (emp-rnf), jpi*jpj    , idex )    ! freshwater budget
721      CALL histwrite( id_i, "sohefldo", kt, qsr + qns, jpi*jpj    , idex )    ! total heat flux
722      CALL histwrite( id_i, "soshfldo", kt, qsr      , jpi*jpj    , idex )    ! solar heat flux
723      CALL histwrite( id_i, "soicecov", kt, fr_i     , jpi*jpj    , idex )    ! ice fraction
724      CALL histwrite( id_i, "sozotaux", kt, utau     , jpi*jpj    , idex )    ! i-wind stress
725      CALL histwrite( id_i, "sometauy", kt, vtau     , jpi*jpj    , idex )    ! j-wind stress
726
727      ! 3. Close the file
728      ! -----------------
729      CALL histclo( id_i )
730#if ! defined key_iomput && ! defined key_dimgout
731      IF( ninist /= 1  ) THEN
732         CALL histclo( nid_T )
733         CALL histclo( nid_U )
734         CALL histclo( nid_V )
735         CALL histclo( nid_W )
736      ENDIF
737#endif
738
739   END SUBROUTINE dia_wri_state
740   !!======================================================================
741END MODULE diawri
Note: See TracBrowser for help on using the repository browser.