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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 2715

Last change on this file since 2715 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk

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