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/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/2011/dev_r2787_LOCEAN3_TRA_TRP/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 2789

Last change on this file since 2789 was 2789, checked in by cetlod, 13 years ago

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