source: trunk/NEMOGCM/NEMO/LIM_SRC_2/limwri_2.F90 @ 3558

Last change on this file since 3558 was 3558, checked in by rblod, 8 years ago

Fix issues when using key_nosignedzeo, see ticket #996

  • Property svn:keywords set to Id
File size: 17.1 KB
Line 
1MODULE limwri_2
2   !!======================================================================
3   !!                     ***  MODULE  limwri_2  ***
4   !!         Ice diagnostics :  write ice output files
5   !!======================================================================
6   !! history :  2.0  ! 2003-08  (C. Ethe)      original code
7   !!            2.0  ! 2004-10  (C. Ethe )     1D configuration
8   !!             -   ! 2009-06  (B. Lemaire )  iom_put + lim_wri_state_2
9   !!-------------------------------------------------------------------
10#if defined key_lim2
11   !!----------------------------------------------------------------------
12   !!   'key_lim2'                                    LIM 2.0 sea-ice model
13   !!----------------------------------------------------------------------
14   !!----------------------------------------------------------------------
15   !!   lim_wri_2      : write of the diagnostics variables in ouput file
16   !!   lim_wri_init_2 : initialization and namelist read
17   !!   lim_wri_state_2 : write for initial state or/and abandon:
18   !!                     > output.init.nc (if ninist = 1 in namelist)
19   !!                     > output.abort.nc
20   !!----------------------------------------------------------------------
21   USE phycst
22   USE dom_oce
23   USE sbc_oce
24   USE sbc_ice
25   USE dom_ice_2
26   USE ice_2
27
28   USE dianam          ! build name of file (routine)
29   USE lbclnk
30   USE in_out_manager
31   USE lib_mpp         ! MPP library
32   USE wrk_nemo        ! work arrays
33   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
34   USE iom
35   USE ioipsl
36
37   IMPLICIT NONE
38   PRIVATE
39
40#if ! defined key_iomput
41   PUBLIC   lim_wri_2         ! called by sbc_ice_lim_2
42#endif
43   PUBLIC   lim_wri_state_2   ! called by dia_wri_state
44   PUBLIC   lim_wri_alloc_2   ! called by nemogcm.F90
45
46   INTEGER, PARAMETER                       ::   jpnoumax = 40   ! maximum number of variable for ice output
47   INTEGER                                  ::   noumef          ! number of fields
48   REAL(wp)           , DIMENSION(jpnoumax) ::   cmulti ,     &  ! multiplicative constant
49      &                                          cadd            ! additive constant
50   CHARACTER(len = 35), DIMENSION(jpnoumax) ::   titn            ! title of the field
51   CHARACTER(len = 8 ), DIMENSION(jpnoumax) ::   nam             ! name of the field
52   CHARACTER(len = 8 ), DIMENSION(jpnoumax) ::   uni             ! unit of the field
53   INTEGER            , DIMENSION(jpnoumax) ::   nc              ! switch for saving field ( = 1 ) or not ( = 0 )
54
55   INTEGER ::   nice, nhorid, ndim, niter, ndepid       ! ????
56   INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndex51   ! ????
57
58   REAL(wp) ::   epsi16 = 1.e-16_wp   ! constant values
59   REAL(wp) ::   zzero  = 0._wp       !     -      -
60   REAL(wp) ::   zone   = 1._wp       !     -      -
61
62   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zcmo      ! Workspace array for netcdf writer.
63
64
65   !! * Substitutions
66#   include "vectopt_loop_substitute.h90"
67   !!----------------------------------------------------------------------
68   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
69   !! $Id$
70   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
71   !!----------------------------------------------------------------------
72CONTAINS
73
74   INTEGER FUNCTION lim_wri_alloc_2()
75      !!-------------------------------------------------------------------
76      !!                  ***   ROUTINE lim_wri_alloc_2  ***
77      !!-------------------------------------------------------------------
78      ALLOCATE( ndex51(jpij), zcmo(jpi,jpj,jpnoumax), STAT=lim_wri_alloc_2)
79      !
80      IF( lk_mpp               )   CALL mpp_sum ( lim_wri_alloc_2 )
81      IF( lim_wri_alloc_2 /= 0 )   CALL ctl_warn('lim_wri_alloc_2: failed to allocate array ndex51')
82      !
83   END FUNCTION lim_wri_alloc_2
84
85
86#if ! defined key_iomput
87# if defined key_dimgout
88   !!----------------------------------------------------------------------
89   !!   'key_dimgout'                                    Direct Access file
90   !!----------------------------------------------------------------------
91# include "limwri_dimg_2.h90"
92# else
93   SUBROUTINE lim_wri_2( kt )
94      !!-------------------------------------------------------------------
95      !!                    ***   ROUTINE lim_wri_2  ***
96      !!               
97      !! ** Purpose :   write the sea-ice output file in NetCDF
98      !!
99      !! ** Method  :   computes the average of some variables and write
100      !!      it in the NetCDF ouput files
101      !!      CAUTION: the sea-ice time-step must be an integer fraction
102      !!      of a day
103      !!-------------------------------------------------------------------
104      INTEGER, INTENT(in) ::   kt     ! number of iteration
105      !!
106      INTEGER  ::   ji, jj, jf                      ! dummy loop indices
107      CHARACTER(len = 40)  ::   clhstnam, clop
108      REAL(wp) ::   zsto, zjulian, zout,   &  ! temporary scalars
109         &          zindh, zinda, zindb, ztmu
110      REAL(wp), DIMENSION(1)                ::   zdept
111      REAL(wp), POINTER, DIMENSION(:,:)     ::   zfield
112      !!-------------------------------------------------------------------
113
114      CALL wrk_alloc( jpi, jpj, zfield )
115                                                 !--------------------!
116      IF( kt == nit000 ) THEN                    !   Initialisation   !
117         !                                       !--------------------!
118
119         CALL lim_wri_init_2 
120                           
121         zsto     = rdt_ice
122         IF( ln_mskland )   THEN   ;   clop = "ave(only(x))"   ! put 1.e+20 on land (very expensive!!)
123         ELSE                      ;   clop = "ave(x)"         ! no use of the mask value (require less cpu time)
124         ENDIF
125         zout     = nwrite * rdt_ice / nn_fsbc
126         niter    = ( nit000 - 1 ) / nn_fsbc
127         zdept(1) = 0.
128         
129         CALL ymds2ju ( nyear, nmonth, nday, rdt, zjulian )
130         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
131         CALL dia_nam ( clhstnam, nwrite, 'icemod' )
132         CALL histbeg ( clhstnam, jpi, glamt, jpj, gphit,    &
133            &           1, jpi, 1, jpj, niter, zjulian, rdt_ice, nhorid, nice , domain_id=nidom, snc4chunks=snc4set)
134         CALL histvert( nice, "deptht", "Vertical T levels", "m", 1, zdept, ndepid, "down")
135         CALL wheneq  ( jpij , tmask(:,:,1), 1, 1., ndex51, ndim)
136         
137         DO jf = 1, noumef
138            IF( nc(jf) == 1 )   CALL histdef( nice, nam(jf), titn(jf), uni(jf), jpi, jpj   &
139               &                                  , nhorid, 1, 1, 1, -99, 32, clop, zsto, zout )
140         END DO
141         CALL histend( nice, snc4set )
142         !
143      ENDIF
144      !                                          !--------------------!
145      !                                          !   Cumulate at kt   !
146      !                                          !--------------------!
147
148      !-- Store instantaneous values in zcmo
149     
150      zcmo(:,:, 1:jpnoumax ) = 0.e0 
151      DO jj = 2 , jpjm1
152         DO ji = 1 , jpim1   ! NO vector opt.
153            zindh  = MAX( zzero , SIGN( zone , hicif(ji,jj) * (1.0 - frld(ji,jj) ) - 0.10 ) )
154            zinda  = MAX( zzero , SIGN( zone , ( 1.0 - frld(ji,jj) ) - 0.10 ) )
155            zindb  = zindh * zinda
156            ztmu   = MAX( 0.5 * zone , ( tmu(ji,jj) + tmu(ji+1,jj) + tmu(ji,jj+1) + tmu(ji+1,jj+1) ) ) 
157            zcmo(ji,jj,1)  = hsnif (ji,jj)
158            zcmo(ji,jj,2)  = hicif (ji,jj)
159            zcmo(ji,jj,3)  = hicifp(ji,jj)
160            zcmo(ji,jj,4)  = frld  (ji,jj)
161            zcmo(ji,jj,5)  = sist  (ji,jj)
162            zcmo(ji,jj,6)  = fbif  (ji,jj)
163            zcmo(ji,jj,7)  = zindb * (  u_ice(ji,jj  ) * tmu(ji,jj  ) + u_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
164                                      + u_ice(ji,jj+1) * tmu(ji,jj+1) + u_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
165                                  / ztmu 
166
167            zcmo(ji,jj,8)  = zindb * (  v_ice(ji,jj  ) * tmu(ji,jj  ) + v_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
168                                      + v_ice(ji,jj+1) * tmu(ji,jj+1) + v_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
169                                  / ztmu
170            zcmo(ji,jj,9)  = sst_m(ji,jj)
171            zcmo(ji,jj,10) = sss_m(ji,jj)
172            zcmo(ji,jj,11) = qns(ji,jj) + qsr(ji,jj)
173            zcmo(ji,jj,12) = qsr(ji,jj)
174            zcmo(ji,jj,13) = qns(ji,jj)
175            ! See thersf for the coefficient
176            zcmo(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce    !!gm ???
177            zcmo(ji,jj,15) = utau_ice(ji,jj)
178            zcmo(ji,jj,16) = vtau_ice(ji,jj)
179            zcmo(ji,jj,17) = qsr_ice(ji,jj,1)
180            zcmo(ji,jj,18) = qns_ice(ji,jj,1)
181            zcmo(ji,jj,19) = sprecip(ji,jj)
182         END DO
183      END DO
184      !
185      ! Write the netcdf file
186      !
187      niter = niter + 1
188      DO jf = 1 , noumef
189         DO jj = 1 , jpj
190            DO ji = 1 , jpi
191               zfield(ji,jj) = zcmo(ji,jj,jf) * cmulti(jf) + cadd(jf)
192            END DO
193         END DO
194         
195         IF( jf == 7  .OR. jf == 8  .OR. jf == 15 .OR. jf == 16 ) THEN
196            CALL lbc_lnk( zfield, 'T', -1. )
197         ELSE
198            CALL lbc_lnk( zfield, 'T',  1. )
199         ENDIF
200         
201         IF( nc(jf) == 1 )   CALL histwrite( nice, nam(jf), niter, zfield, ndim, ndex51 )
202         
203      END DO
204     
205      IF( ( nn_fsbc * niter ) >= nitend )   CALL histclo( nice ) 
206
207      CALL wrk_dealloc( jpi, jpj, zfield )
208      !
209   END SUBROUTINE lim_wri_2
210     
211# endif
212
213   SUBROUTINE lim_wri_init_2
214      !!-------------------------------------------------------------------
215      !!                    ***   ROUTINE lim_wri_init_2  ***
216      !!               
217      !! ** Purpose :   intialisation of LIM sea-ice output
218      !!
219      !! ** Method  : Read the namicewri namelist and check the parameter
220      !!       values called at the first timestep (nit000)
221      !!
222      !! ** input   :   Namelist namicewri
223      !!-------------------------------------------------------------------
224      INTEGER ::   nf      ! ???
225      TYPE FIELD 
226         CHARACTER(len = 35) :: ztitle 
227         CHARACTER(len = 8 ) :: zname         
228         CHARACTER(len = 8 ) :: zunit
229         INTEGER             :: znc   
230         REAL                :: zcmulti 
231         REAL                :: zcadd       
232      END TYPE FIELD
233      TYPE(FIELD) ::  &
234         field_1 , field_2 , field_3 , field_4 , field_5 , field_6 ,   &
235         field_7 , field_8 , field_9 , field_10, field_11, field_12,   &
236         field_13, field_14, field_15, field_16, field_17, field_18,   &
237         field_19
238      TYPE(FIELD) , DIMENSION(jpnoumax) :: zfield
239
240      NAMELIST/namiceout/ noumef, &
241         field_1 , field_2 , field_3 , field_4 , field_5 , field_6 ,   &
242         field_7 , field_8 , field_9 , field_10, field_11, field_12,   &
243         field_13, field_14, field_15, field_16, field_17, field_18,   &
244         field_19
245      !!-------------------------------------------------------------------
246      !
247      IF( lim_wri_alloc_2() /= 0 ) THEN      ! allocate lim_wri arrrays
248         CALL ctl_stop( 'STOP', 'lim_wri_init_2 : unable to allocate standard arrays' )   ;   RETURN
249      ENDIF
250
251      REWIND ( numnam_ice )                ! Read Namelist namicewri
252      READ   ( numnam_ice  , namiceout )
253     
254      zfield( 1) = field_1
255      zfield( 2) = field_2
256      zfield( 3) = field_3
257      zfield( 4) = field_4
258      zfield( 5) = field_5
259      zfield( 6) = field_6
260      zfield( 7) = field_7
261      zfield( 8) = field_8
262      zfield( 9) = field_9
263      zfield(10) = field_10
264      zfield(11) = field_11
265      zfield(12) = field_12
266      zfield(13) = field_13
267      zfield(14) = field_14
268      zfield(15) = field_15
269      zfield(16) = field_16
270      zfield(17) = field_17
271      zfield(18) = field_18
272      zfield(19) = field_19
273     
274      DO nf = 1, noumef
275         titn  (nf) = zfield(nf)%ztitle
276         nam   (nf) = zfield(nf)%zname
277         uni   (nf) = zfield(nf)%zunit
278         nc    (nf) = zfield(nf)%znc
279         cmulti(nf) = zfield(nf)%zcmulti
280         cadd  (nf) = zfield(nf)%zcadd
281      END DO
282
283      IF(lwp) THEN
284         WRITE(numout,*)
285         WRITE(numout,*) 'lim_wri_init_2 : Ice parameters for outputs'
286         WRITE(numout,*) '~~~~~~~~~~~~~~'
287         WRITE(numout,*) '    number of fields to be stored         noumef = ', noumef
288         WRITE(numout,*) '           title                            name     unit      Saving (1/0) ',   &
289            &            '    multiplicative constant       additive constant '
290         DO nf = 1 , noumef         
291            WRITE(numout,*) '   ', titn(nf), '   ', nam(nf),'      ', uni(nf),'  ', nc(nf),'        ', cmulti(nf),   &
292               &       '        ', cadd(nf)
293         END DO
294      ENDIF
295      !   
296   END SUBROUTINE lim_wri_init_2
297
298#endif
299
300   SUBROUTINE lim_wri_state_2( kt, kid, kh_i )
301      !!---------------------------------------------------------------------
302      !!                 ***  ROUTINE lim_wri_state_2  ***
303      !!       
304      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
305      !!      the instantaneous ice state and forcing fields for ice model
306      !!        Used to find errors in the initial state or save the last
307      !!      ocean state in case of abnormal end of a simulation
308      !!
309      !! History :
310      !!   2.0  !  2009-06  (B. Lemaire)
311      !!----------------------------------------------------------------------
312      INTEGER, INTENT( in ) ::   kt               ! ocean time-step index)
313      INTEGER, INTENT( in ) ::   kid , kh_i       
314      !!----------------------------------------------------------------------
315
316      CALL histdef( kid, "isnowthi", "Snow thickness"          , "m"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
317      CALL histdef( kid, "iicethic", "Ice thickness"           , "m"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
318      CALL histdef( kid, "iiceprod", "Ice produced"            , "m/kt"   , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
319      CALL histdef( kid, "ileadfra", "Ice concentration"       , "-"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
320      CALL histdef( kid, "iicetemp", "Ice temperature"         , "K"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
321      CALL histdef( kid, "ioceflxb", "flux at ice base"        , "w/m2"   , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
322      CALL histdef( kid, "iicevelu", "i-Ice speed (I-point)"   , "m/s"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
323      CALL histdef( kid, "iicevelv", "j-Ice speed (I-point)"   , "m/s"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
324      CALL histdef( kid, "isstempe", "Sea surface temperature" , "C"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
325      CALL histdef( kid, "isssalin", "Sea surface salinity"    , "PSU"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
326      CALL histdef( kid, "iicestru", "i-Wind stress over ice (I-pt)", "Pa", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
327      CALL histdef( kid, "iicestrv", "j-Wind stress over ice (I-pt)", "Pa", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
328      CALL histdef( kid, "iicesflx", "Solar flux over ice"     , "w/m2"   , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
329      CALL histdef( kid, "iicenflx", "Non-solar flux over ice" , "w/m2"   , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
330      CALL histdef( kid, "isnowpre", "Snow precipitation"      , "kg/m2/s", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
331
332      CALL histend( kid, snc4set )   ! end of the file definition
333
334      CALL histwrite( kid, "isnowthi", kt, hsnif          , jpi*jpj, (/1/) )   
335      CALL histwrite( kid, "iicethic", kt, hicif          , jpi*jpj, (/1/) )   
336      CALL histwrite( kid, "iiceprod", kt, hicifp         , jpi*jpj, (/1/) )   
337      CALL histwrite( kid, "ileadfra", kt, 1. - frld(:,:) , jpi*jpj, (/1/) )
338      CALL histwrite( kid, "iicetemp", kt, sist(:,:) - rt0, jpi*jpj, (/1/) )
339      CALL histwrite( kid, "ioceflxb", kt, fbif           , jpi*jpj, (/1/) )
340      CALL histwrite( kid, "iicevelu", kt, u_ice          , jpi*jpj, (/1/) )
341      CALL histwrite( kid, "iicevelv", kt, v_ice          , jpi*jpj, (/1/) )
342      CALL histwrite( kid, "isstempe", kt, sst_m          , jpi*jpj, (/1/) )
343      CALL histwrite( kid, "isssalin", kt, sss_m          , jpi*jpj, (/1/) )
344      CALL histwrite( kid, "iicestru", kt, utau_ice       , jpi*jpj, (/1/) )
345      CALL histwrite( kid, "iicestrv", kt, vtau_ice       , jpi*jpj, (/1/) )
346      CALL histwrite( kid, "iicesflx", kt, qsr_ice(:,:,1) , jpi*jpj, (/1/) )
347      CALL histwrite( kid, "iicenflx", kt, qns_ice(:,:,1) , jpi*jpj, (/1/) )
348      CALL histwrite( kid, "isnowpre", kt, sprecip        , jpi*jpj, (/1/) )
349
350    END SUBROUTINE lim_wri_state_2
351
352#else
353   !!----------------------------------------------------------------------
354   !!   Default option :         Empty module      NO LIM 2.0 sea-ice model
355   !!----------------------------------------------------------------------
356CONTAINS
357   SUBROUTINE lim_wri_2          ! Empty routine
358   END SUBROUTINE lim_wri_2
359#endif
360
361   !!======================================================================
362END MODULE limwri_2
Note: See TracBrowser for help on using the repository browser.