source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90 @ 8413

Last change on this file since 8413 was 8413, checked in by clem, 3 years ago

continue changing names (again)

File size: 36.5 KB
Line 
1MODULE icewri
2   !!======================================================================
3   !!                     ***  MODULE  icewri  ***
4   !!         Ice diagnostics :  write ice output files
5   !!======================================================================
6#if defined key_lim3
7   !!----------------------------------------------------------------------
8   !!   'key_lim3'                                      LIM3 sea-ice model
9   !!----------------------------------------------------------------------
10   !!   ice_wri      : write of the diagnostics variables in ouput file
11   !!   ice_wri_state : write for initial state or/and abandon
12   !!----------------------------------------------------------------------
13   USE ioipsl
14   USE dianam          ! build name of file (routine)
15   USE phycst
16   USE dom_oce
17   USE sbc_oce
18   USE sbc_ice         ! Surface boundary condition: ice fields
19   USE ice
20   USE limvar
21   !
22   USE in_out_manager
23   USE lbclnk
24   USE lib_mpp         ! MPP library
25   USE iom
26   USE timing          ! Timing
27   USE lib_fortran     ! Fortran utilities
28
29   IMPLICIT NONE
30   PRIVATE
31
32   PUBLIC ice_wri        ! routine called by lim_step.F90
33   PUBLIC ice_wri_state  ! called by dia_wri_state
34
35   !!----------------------------------------------------------------------
36   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
37   !! $Id: icewri.F90 8409 2017-08-07 15:29:21Z clem $
38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE ice_wri( kindic )
43      !!-------------------------------------------------------------------
44      !!  This routine computes the average of some variables and write it
45      !!  on the ouput files.
46      !!  ATTENTION cette routine n'est valable que si le pas de temps est
47      !!  egale a une fraction entiere de 1 jours.
48      !!  Diff 1-D 3-D : suppress common also included in etat
49      !!                 suppress cmoymo 11-18
50      !!  modif : 03/06/98
51      !!-------------------------------------------------------------------
52      INTEGER, INTENT(in) ::   kindic   ! if kindic < 0 there has been an error somewhere
53      !
54      INTEGER  ::  ji, jj, jk, jl  ! dummy loop indices
55      REAL(wp) ::  z2da, z2db, ztmp, zrho1, zrho2, zmiss_val
56      REAL(wp) ::  zs12, zshear
57      REAL(wp), DIMENSION(jpi,jpj,jpl) ::  zswi2, zmiss2
58      REAL(wp), DIMENSION(jpi,jpj)   ::  z2d, zswi, zmiss !  2D workspace
59      REAL(wp), DIMENSION(jpi,jpj)   ::  zfb              ! ice freeboard
60      REAL(wp), DIMENSION(jpi,jpj)   ::  zamask, zamask15 ! 15% concentration mask
61      REAL(wp), DIMENSION(jpi,jpj)   ::  zsig1, zsig2, zsig3
62
63      ! Global ice diagnostics (SIMIP)
64      REAL(wp) ::  zdiag_area_nh, &   ! area, extent, volume
65         &         zdiag_extt_nh, &
66         &         zdiag_area_sh, & 
67         &         zdiag_extt_sh, & 
68         &         zdiag_volu_nh, & 
69         &         zdiag_volu_sh 
70
71      !!-------------------------------------------------------------------
72
73      IF( nn_timing == 1 )  CALL timing_start('icewri')
74
75      !----------------------------------------
76      ! Brine volume, switches, missing values
77      !----------------------------------------
78
79      ! brine volume
80      CALL lim_var_bv 
81
82      ! tresholds for outputs
83      DO jj = 1, jpj
84         DO ji = 1, jpi
85            zswi(ji,jj)      = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
86            zamask(ji,jj)    = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.05   ) ) ! 1 if 5% ice, 0 if less - required to mask thickness and snow depth
87            zamask15(ji,jj)  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15   ) ) ! 1 if 15% ice, 0 if less
88         END DO
89      END DO
90      DO jl = 1, jpl
91         DO jj = 1, jpj
92            DO ji = 1, jpi
93               zswi2(ji,jj,jl)  = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) )
94            END DO
95         END DO
96      END DO
97
98      zmiss_val     = 1.0e20
99      zmiss(:,:)    = zmiss_val * ( 1. - zswi(:,:) )
100      zmiss2(:,:,:) = zmiss_val * ( 1. - zswi2(:,:,:) )
101
102      !----------------------------------------
103      ! Standard outputs
104      !----------------------------------------
105      ! fluxes
106      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * ( 1._wp - at_i_b(:,:) ) )                      !     solar flux at ocean surface
107      IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + qemp_oce(:,:) )      ! non-solar flux at ocean surface
108      IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface
109      IF( iom_use('qns_ice') )   CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface
110      IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice
111      IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * ( 1._wp - at_i_b(:,:) ) + qemp_oce(:,:) )
112      IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   &
113         &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) )
114      IF( iom_use('qemp_oce') )  CALL iom_put( "qemp_oce" , qemp_oce(:,:) ) 
115      IF( iom_use('qemp_ice') )  CALL iom_put( "qemp_ice" , qemp_ice(:,:) ) 
116      IF( iom_use('emp_oce' ) )  CALL iom_put( "emp_oce"  , emp_oce(:,:) )   ! emp over ocean (taking into account the snow blown away from the ice)
117      IF( iom_use('emp_ice' ) )  CALL iom_put( "emp_ice"  , emp_ice(:,:) )   ! emp over ice   (taking into account the snow blown away from the ice)
118
119      ! velocity
120      IF ( iom_use('uice_ipa') ) CALL iom_put( "uice_ipa"     , u_ice      )       ! ice velocity u component
121      IF ( iom_use('vice_ipa') ) CALL iom_put( "vice_ipa"     , v_ice      )       ! ice velocity v component
122
123      IF ( ( iom_use( "icevel" ) ) .OR. ( iom_use( "icevel_mv" ) ) ) THEN
124         DO jj = 2 , jpjm1
125            DO ji = 2 , jpim1
126               z2da  = ( u_ice(ji,jj) * umask(ji,jj,1) + u_ice(ji-1,jj) * umask(ji-1,jj,1) ) * 0.5_wp
127               z2db  = ( v_ice(ji,jj) * vmask(ji,jj,1) + v_ice(ji,jj-1) * vmask(ji,jj-1,1) ) * 0.5_wp
128               z2d(ji,jj) = SQRT( z2da * z2da + z2db * z2db )
129           END DO
130         END DO
131         CALL lbc_lnk( z2d, 'T', 1. )
132         IF ( iom_use( "icevel" )  )   CALL iom_put( "icevel"       , z2d        )                                            ! ice velocity module
133         IF ( iom_use( "icevel_mv" ) ) CALL iom_put( "icevel_mv"    , z2d(:,:) * zswi(:,:) + zmiss(:,:) )                     ! ice velocity module (missing value)
134      ENDIF
135
136      IF ( iom_use( "tau_icebfr" ) )    CALL iom_put( "tau_icebfr"  , tau_icebfr             )  ! ice friction with ocean bottom (landfast ice) 
137      !
138      IF ( iom_use( "miceage" ) )       CALL iom_put( "miceage"     , om_i * zswi * zamask15 )  ! mean ice age
139      IF ( iom_use( "micet" ) )         CALL iom_put( "micet"       , ( tm_i  - rt0 ) * zswi )  ! ice mean    temperature
140      IF ( iom_use( "icest" ) )         CALL iom_put( "icest"       , ( tm_su - rt0 ) * zswi )  ! ice surface temperature
141      IF ( iom_use( "icecolf" ) )       CALL iom_put( "icecolf"     , hicol                  )  ! frazil ice collection thickness
142      !
143      CALL iom_put( "isst"        , sst_m               )        ! sea surface temperature
144      CALL iom_put( "isss"        , sss_m               )        ! sea surface salinity
145      CALL iom_put( "iceconc"     , at_i  * zswi        )        ! ice concentration
146      CALL iom_put( "icevolu"     , vt_i  * zswi        )        ! ice volume = mean ice thickness over the cell
147      CALL iom_put( "icethick"    , htm_i * zswi        )        ! ice thickness
148      CALL iom_put( "icehc"       , et_i  * zswi        )        ! ice total heat content
149      CALL iom_put( "isnowhc"     , et_s  * zswi        )        ! snow total heat content
150      CALL iom_put( "ibrinv"      , bvm_i * zswi * 100. )        ! brine volume
151      CALL iom_put( "snowpre"     , sprecip * 86400.    )        ! snow precipitation
152      CALL iom_put( "micesalt"    , smt_i   * zswi      )        ! mean ice salinity
153      CALL iom_put( "snowvol"     , vt_s    * zswi      )        ! snow volume
154
155      CALL iom_put( "idive"       , divu_i(:,:)  * zswi(:,:) + zmiss(:,:) )   ! divergence
156      CALL iom_put( "ishear"      , shear_i(:,:) * zswi(:,:) + zmiss(:,:) )   ! shear
157     
158      CALL iom_put( "icetrp"      , diag_trp_vi * rday  )        ! ice volume transport
159      CALL iom_put( "snwtrp"      , diag_trp_vs * rday  )        ! snw volume transport
160      CALL iom_put( "saltrp"      , diag_trp_smv * rday * rhoic ) ! salt content transport
161      CALL iom_put( "deitrp"      , diag_trp_ei         )        ! advected ice enthalpy (W/m2)
162      CALL iom_put( "destrp"      , diag_trp_es         )        ! advected snw enthalpy (W/m2)
163
164      CALL iom_put( "sfxbog"      , sfx_bog * rday      )        ! salt flux from bottom growth
165      CALL iom_put( "sfxbom"      , sfx_bom * rday      )        ! salt flux from bottom melting
166      CALL iom_put( "sfxsum"      , sfx_sum * rday      )        ! salt flux from surface melting
167      CALL iom_put( "sfxlam"      , sfx_lam * rday      )        ! salt flux from lateral melting
168      CALL iom_put( "sfxsni"      , sfx_sni * rday      )        ! salt flux from snow ice formation
169      CALL iom_put( "sfxopw"      , sfx_opw * rday      )        ! salt flux from open water formation
170      CALL iom_put( "sfxdyn"      , sfx_dyn * rday      )        ! salt flux from ridging rafting
171      CALL iom_put( "sfxres"      , sfx_res * rday      )        ! salt flux from corrections (resultant)
172      CALL iom_put( "sfxbri"      , sfx_bri * rday      )        ! salt flux from brines
173      CALL iom_put( "sfxsub"      , sfx_sub * rday      )        ! salt flux from sublimation
174      CALL iom_put( "sfx"         , sfx     * rday      )        ! total salt flux
175
176      ztmp = rday / rhoic
177      CALL iom_put( "vfxres"     , wfx_res * ztmp       )        ! daily prod./melting due to corrections
178      CALL iom_put( "vfxopw"     , wfx_opw * ztmp       )        ! daily lateral thermodynamic ice production
179      CALL iom_put( "vfxsni"     , wfx_sni * ztmp       )        ! daily snowice ice production
180      CALL iom_put( "vfxbog"     , wfx_bog * ztmp       )        ! daily bottom thermodynamic ice production
181      CALL iom_put( "vfxdyn"     , wfx_dyn * ztmp       )        ! daily dynamic ice production (rid/raft)
182      CALL iom_put( "vfxsum"     , wfx_sum * ztmp       )        ! surface melt
183      CALL iom_put( "vfxbom"     , wfx_bom * ztmp       )        ! bottom melt
184      CALL iom_put( "vfxlam"     , wfx_lam * ztmp       )        ! lateral melt
185      CALL iom_put( "vfxice"     , wfx_ice * ztmp       )        ! total ice growth/melt
186
187      IF ( ln_pnd ) &
188         CALL iom_put( "vfxpnd"  , wfx_pnd * ztmp       )        ! melt pond water flux
189
190      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations 
191         WHERE( htm_i(:,:) < 0.2 .AND. htm_i(:,:) > 0. ) ; z2d = wfx_bog
192         ELSEWHERE                                       ; z2d = 0._wp
193         END WHERE
194         CALL iom_put( "vfxthin", ( wfx_opw + z2d ) * ztmp )
195      ENDIF
196
197      ztmp = rday / rhosn
198      CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow)
199      CALL iom_put( "vfxsnw"     , wfx_snw * ztmp       )        ! total snw growth/melt
200      CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow/ice)
201      CALL iom_put( "vfxsub_err" , wfx_err_sub * ztmp   )        ! "excess" of sublimation sent to ocean     
202 
203      CALL iom_put( "afxtot"     , afx_tot              )        ! concentration tendency (total)
204      CALL iom_put( "afxdyn"     , afx_dyn              )        ! concentration tendency (dynamics)
205      CALL iom_put( "afxthd"     , afx_thd              )        ! concentration tendency (thermo)
206
207      CALL iom_put ('hfxthd'     , hfx_thd(:,:)         )   
208      CALL iom_put ('hfxdyn'     , hfx_dyn(:,:)         )   
209      CALL iom_put ('hfxres'     , hfx_res(:,:)         )   
210      CALL iom_put ('hfxout'     , hfx_out(:,:)         )   
211      CALL iom_put ('hfxin'      , hfx_in(:,:)          )   
212      CALL iom_put ('hfxsnw'     , hfx_snw(:,:)         )   
213      CALL iom_put ('hfxsub'     , hfx_sub(:,:)         )   
214      CALL iom_put ('hfxerr'     , hfx_err(:,:)         )   
215      CALL iom_put ('hfxerr_rem' , hfx_err_rem(:,:)     )   
216     
217      CALL iom_put ('hfxsum'     , hfx_sum(:,:)         )   
218      CALL iom_put ('hfxbom'     , hfx_bom(:,:)         )   
219      CALL iom_put ('hfxbog'     , hfx_bog(:,:)         )   
220      CALL iom_put ('hfxdif'     , hfx_dif(:,:)         )   
221      CALL iom_put ('hfxopw'     , hfx_opw(:,:)         )   
222      CALL iom_put ('hfxtur'     , fhtur(:,:) * at_i_b(:,:) ) ! turbulent heat flux at ice base
223      CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice
224      CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip
225
226      ! specific outputs for EVP rheology
227      IF( iom_use( "isig1" ) .OR. iom_use( "isig2" ) .OR. iom_use( "isig3" ) ) THEN
228         zsig1(:,:) = 0._wp; zsig2(:,:) = 0._wp; zsig3(:,:) = 0._wp;
229         DO jj = 2, jpjm1
230            DO ji = 2, jpim1
231               zs12 = ( zswi(ji-1,jj) * stress12_i(ji-1,jj) + zswi(ji  ,jj-1) * stress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
232                  &     zswi(ji  ,jj) * stress12_i(ji  ,jj) + zswi(ji-1,jj-1) * stress12_i(ji-1,jj-1) )  &
233                  &   / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) )
234
235               zshear = SQRT( stress2_i(ji,jj) * stress2_i(ji,jj) + 4._wp * zs12 * zs12 ) ! shear stress 
236
237               z2da = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) )
238
239!!               zsig1(ji,jj) = 0.5_wp * z2da * ( stress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
240!!               zsig2(ji,jj) = 0.5_wp * z2da * ( stress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
241!!               zsig3(ji,jj) = z2da**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
242!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
243               zsig1(ji,jj) = 0.5_wp * z2da * ( stress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
244               zsig2(ji,jj) = 0.5_wp * z2da * ( zshear )                    ! shear stress
245               zsig3(ji,jj) = z2da**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
246            END DO
247         END DO
248         CALL lbc_lnk_multi( zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
249         CALL iom_put( "isig1" , zsig1 )
250         CALL iom_put( "isig2" , zsig2 )
251         CALL iom_put( "isig3" , zsig3 )
252      ENDIF
253
254      ! MV MP 2016
255      IF ( ln_pnd ) THEN
256         CALL iom_put( "iceamp"  , at_ip  * zswi        )   ! melt pond total fraction
257         CALL iom_put( "icevmp"  , vt_ip  * zswi        )   ! melt pond total volume per unit area
258      ENDIF
259      ! END MV MP 2016
260
261      !----------------------------------
262      ! Output category-dependent fields
263      !----------------------------------
264      IF ( iom_use( "iceconc_cat"  ) )  CALL iom_put( "iceconc_cat"      , a_i   * zswi2   )        ! area for categories
265      IF ( iom_use( "icethic_cat"  ) )  CALL iom_put( "icethic_cat"      , ht_i  * zswi2   )        ! thickness for categories
266      IF ( iom_use( "snowthic_cat" ) )  CALL iom_put( "snowthic_cat"     , ht_s  * zswi2   )        ! snow depth for categories
267      IF ( iom_use( "salinity_cat" ) )  CALL iom_put( "salinity_cat"     , sm_i  * zswi2   )        ! salinity for categories
268      ! ice temperature
269      IF ( iom_use( "icetemp_cat"  ) )  CALL iom_put( "icetemp_cat", ( SUM( t_i(:,:,:,:), dim=3 ) * r1_nlay_i - rt0 ) * zswi2 )
270      ! snow temperature
271      IF ( iom_use( "snwtemp_cat"  ) )  CALL iom_put( "snwtemp_cat", ( SUM( t_s(:,:,:,:), dim=3 ) * r1_nlay_s - rt0 ) * zswi2 )
272      ! ice age
273      IF ( iom_use( "iceage_cat"   ) )  CALL iom_put( "iceage_cat" , o_i * zswi2 ) 
274      ! brine volume
275      IF ( iom_use( "brinevol_cat" ) )  CALL iom_put( "brinevol_cat", bv_i * 100. * zswi2 )
276
277      ! MV MP 2016
278      IF ( ln_pnd ) THEN
279         IF ( iom_use( "iceamp_cat"  ) )  CALL iom_put( "iceamp_cat"     , a_ip       * zswi2   )       ! melt pond frac for categories
280         IF ( iom_use( "icevmp_cat"  ) )  CALL iom_put( "icevmp_cat"     , v_ip       * zswi2   )       ! melt pond frac for categories
281         IF ( iom_use( "icehmp_cat"  ) )  CALL iom_put( "icehmp_cat"     , h_ip       * zswi2   )       ! melt pond frac for categories
282         IF ( iom_use( "iceafp_cat"  ) )  CALL iom_put( "iceafp_cat"     , a_ip_frac  * zswi2   )       ! melt pond frac for categories
283      ENDIF
284      ! END MV MP 2016
285
286      !--------------------------------
287      ! Add-ons for SIMIP
288      !--------------------------------
289      zrho1 = ( rau0 - rhoic ) / rau0; zrho2 = rhosn / rau0
290
291      IF  ( iom_use( "icepres"  ) ) CALL iom_put( "icepres"     , zswi(:,:)                     )                                ! Ice presence (1 or 0)
292      IF  ( iom_use( "icemass"  ) ) CALL iom_put( "icemass"     , rhoic * vt_i(:,:) * zswi(:,:) )                                ! Ice mass per cell area
293      IF  ( iom_use( "icethic"  ) ) CALL iom_put( "icethic"     , htm_i(:,:) * zamask(:,:)  + ( 1. - zamask(:,:) ) * zmiss_val )     ! Ice thickness
294      IF  ( iom_use( "snomass"  ) ) CALL iom_put( "snomass"     , rhosn * vt_s(:,:)         * zswi(:,:) + zmiss(:,:) )           ! Snow mass per cell area
295      IF  ( iom_use( "snothic"  ) ) CALL iom_put( "snothic"     , htm_s(:,:) * zamask(:,:)  + ( 1. - zamask(:,:) ) * zmiss_val )     ! Snow thickness       
296
297      IF  ( iom_use( "iceconc_cat_mv"  ) )  CALL iom_put( "iceconc_cat_mv" , a_i(:,:,:)  * zswi2(:,:,:) + zmiss2(:,:,:) )        ! Area for categories
298      IF  ( iom_use( "icethic_cat_mv"  ) )  CALL iom_put( "icethic_cat_mv" , ht_i(:,:,:) * zswi2(:,:,:) + zmiss2(:,:,:) )        ! Thickness for categories
299      IF  ( iom_use( "snowthic_cat_mv" ) )  CALL iom_put( "snowthic_cat_mv", ht_s(:,:,:) * zswi2(:,:,:) + zmiss2(:,:,:) )        ! Snow depth for categories
300
301      IF  ( iom_use( "icestK"   ) ) CALL iom_put( "icestK"      , tm_su(:,:)                * zswi(:,:) + zmiss(:,:) )           ! Ice surface temperature
302      IF  ( iom_use( "icesntK"  ) ) CALL iom_put( "icesntK"     , tm_si(:,:)                * zswi(:,:) + zmiss(:,:) )           ! Snow-ice interface temperature
303      IF  ( iom_use( "icebotK"  ) ) CALL iom_put( "icebotK"     , t_bo(:,:)                 * zswi(:,:) + zmiss(:,:) )           ! Ice bottom temperature
304      IF  ( iom_use( "iceage"   ) ) CALL iom_put( "iceage"      , om_i(:,:) * zamask15(:,:) + ( 1. - zamask15(:,:) ) * zmiss_val )   ! Ice age
305      IF  ( iom_use( "icesmass" ) ) CALL iom_put( "icesmass"    , SUM( smv_i, DIM = 3 ) * rhoic * 1.0e-3 * zswi(:,:) )           ! Mass of salt in sea ice per cell area
306      IF  ( iom_use( "icesal"   ) ) CALL iom_put( "icesal"      , smt_i(:,:)                * zswi(:,:) + zmiss(:,:) )           ! Ice salinity
307
308      IF  ( iom_use( "icefb"    ) ) THEN
309         zfb(:,:) = ( zrho1 * htm_i(:,:) - zrho2 * htm_s(:,:) )                                         
310         WHERE( zfb < 0._wp ) ;   zfb = 0._wp ;   END WHERE
311                                    CALL iom_put( "icefb"       , zfb(:,:)                  * zswi(:,:) + zmiss(:,:) )           ! Ice freeboard
312      ENDIF
313
314      IF  ( iom_use( "isnhcneg" ) ) CALL iom_put( "isnhcneg"    , - et_s(:,:)               * zswi(:,:) + zmiss(:,:) )           ! Snow total heat content
315
316      IF  ( iom_use( "dmithd"   ) ) CALL iom_put( "dmithd"      , - wfx_bog - wfx_bom - wfx_sum   &                       ! Sea-ice mass change from thermodynamics
317              &                     - wfx_sni - wfx_opw - wfx_res )
318      IF  ( iom_use( "dmidyn"   ) ) CALL iom_put( "dmidyn"      ,   diag_dmi_dyn             )                            ! Sea-ice mass change from dynamics
319      IF  ( iom_use( "dmiopw"   ) ) CALL iom_put( "dmiopw"      , - wfx_opw                  )                            ! Sea-ice mass change through growth in open water
320      IF  ( iom_use( "dmibog"   ) ) CALL iom_put( "dmibog"      , - wfx_bog                  )                            ! Sea-ice mass change through basal growth
321      IF  ( iom_use( "dmisni"   ) ) CALL iom_put( "dmisni"      , - wfx_sni                  )                            ! Sea-ice mass change through snow-to-ice conversion
322      IF  ( iom_use( "dmisum"   ) ) CALL iom_put( "dmisum"      , - wfx_sum                  )                            ! Sea-ice mass change through surface melting
323      IF  ( iom_use( "dmibom"   ) ) CALL iom_put( "dmibom"      , - wfx_bom                  )                            ! Sea-ice mass change through bottom melting
324
325      IF  ( iom_use( "dmtsub"   ) ) CALL iom_put( "dmtsub"      , - wfx_sub                  )                            ! Sea-ice mass change through evaporation and sublimation
326      IF  ( iom_use( "dmssub"   ) ) CALL iom_put( "dmssub"      , - wfx_snw_sub              )                            ! Snow mass change through sublimation
327      IF  ( iom_use( "dmisub"   ) ) CALL iom_put( "dmisub"      , - wfx_ice_sub              )                            ! Sea-ice mass change through sublimation
328
329      IF  ( iom_use( "dmsspr"   ) ) CALL iom_put( "dmsspr"      , - wfx_spr                  )                            ! Snow mass change through snow fall
330      IF  ( iom_use( "dmsssi"   ) ) CALL iom_put( "dmsssi"      ,   wfx_sni*rhosn/rhoic      )                            ! Snow mass change through snow-to-ice conversion
331
332      IF  ( iom_use( "dmsmel"   ) ) CALL iom_put( "dmsmel"      , - wfx_snw_sum              )                            ! Snow mass change through melt
333      IF  ( iom_use( "dmsdyn"   ) ) CALL iom_put( "dmsdyn"      ,   diag_dms_dyn             )                            ! Snow mass change through dynamics
334
335      IF  ( iom_use( "hfxsenso" ) ) CALL iom_put( "hfxsenso"    ,   -fhtur(:,:)              * zswi(:,:) + zmiss(:,:) )   ! Sensible oceanic heat flux
336      IF  ( iom_use( "hfxconbo" ) ) CALL iom_put( "hfxconbo"    ,   diag_fc_bo               * zswi(:,:) + zmiss(:,:) )   ! Bottom conduction flux
337      IF  ( iom_use( "hfxconsu" ) ) CALL iom_put( "hfxconsu"    ,   diag_fc_su               * zswi(:,:) + zmiss(:,:) )   ! Surface conduction flux
338
339      IF  ( iom_use( "wfxtot"   ) ) CALL iom_put( "wfxtot"      ,   wfx_ice(:,:)             * zswi(:,:) + zmiss(:,:) )   ! Total freshwater flux from sea ice
340      IF  ( iom_use( "wfxsum"   ) ) CALL iom_put( "wfxsum"      ,   wfx_sum(:,:)             * zswi(:,:) + zmiss(:,:) )   ! Freshwater flux from sea-ice surface
341      IF  ( iom_use( "sfx_mv"   ) ) CALL iom_put( "sfx_mv"      ,   sfx(:,:) * 0.001         * zswi(:,:) + zmiss(:,:) )   ! Total salt flux
342
343      IF  ( iom_use( "uice_mv"  ) ) CALL iom_put( "uice_mv"     ,   u_ice(:,:)               * zswi(:,:) + zmiss(:,:) )   ! ice velocity u component
344      IF  ( iom_use( "vice_mv"  ) ) CALL iom_put( "vice_mv"     ,   v_ice(:,:)               * zswi(:,:) + zmiss(:,:) )   ! ice velocity v component
345     
346      IF  ( iom_use( "xmtrpice" ) ) CALL iom_put( "xmtrpice"     ,  diag_xmtrp_ice(:,:)      )                            ! X-component of sea-ice mass transport (kg/s)
347      IF  ( iom_use( "ymtrpice" ) ) CALL iom_put( "ymtrpice"     ,  diag_ymtrp_ice(:,:)      )                            ! Y-component of sea-ice mass transport
348
349      IF  ( iom_use( "xmtrpsnw" ) ) CALL iom_put( "xmtrpsnw"     ,  diag_xmtrp_snw(:,:)      )                            ! X-component of snow mass transport (kg/s)
350      IF  ( iom_use( "ymtrpsnw" ) ) CALL iom_put( "ymtrpsnw"     ,  diag_ymtrp_snw(:,:)      )                            ! Y-component of snow mass transport
351
352      IF  ( iom_use( "xatrp"    ) ) CALL iom_put( "xatrp"        ,  diag_xatrp(:,:)              )                        ! X-component of ice area transport
353      IF  ( iom_use( "yatrp"    ) ) CALL iom_put( "yatrp"        ,  diag_yatrp(:,:)              )                        ! Y-component of ice area transport
354
355      IF  ( iom_use( "utau_ice" ) ) CALL iom_put( "utau_ice"     ,  utau_ice(:,:)            * zswi(:,:) + zmiss(:,:) )   ! Wind stress term in force balance (x)
356      IF  ( iom_use( "vtau_ice" ) ) CALL iom_put( "vtau_ice"     ,  vtau_ice(:,:)            * zswi(:,:) + zmiss(:,:) )   ! Wind stress term in force balance (y)
357
358      IF  ( iom_use( "utau_oi"  ) ) CALL iom_put( "utau_oi"     ,   diag_utau_oi(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Ocean stress term in force balance (x)
359      IF  ( iom_use( "vtau_oi"  ) ) CALL iom_put( "vtau_oi"     ,   diag_vtau_oi(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Ocean stress term in force balance (y)
360
361      IF  ( iom_use( "icestr"   ) ) CALL iom_put( "icestr"      ,   strength(:,:)            * zswi(:,:) + zmiss(:,:) )   ! Ice strength
362
363      IF  ( iom_use( "dssh_dx"  ) ) CALL iom_put( "dssh_dx"     ,   diag_dssh_dx(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Sea-surface tilt term in force balance (x)
364      IF  ( iom_use( "dssh_dy"  ) ) CALL iom_put( "dssh_dy"     ,   diag_dssh_dy(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Sea-surface tilt term in force balance (y)
365
366      IF  ( iom_use( "corstrx"  ) ) CALL iom_put( "corstrx"     ,   diag_corstrx(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Coriolis force term in force balance (x)
367      IF  ( iom_use( "corstry"  ) ) CALL iom_put( "corstry"     ,   diag_corstry(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Coriolis force term in force balance (y)
368
369      IF  ( iom_use( "intstrx"  ) ) CALL iom_put( "intstrx"     ,   diag_intstrx(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Internal force term in force balance (x)
370      IF  ( iom_use( "intstry"  ) ) CALL iom_put( "intstry"     ,   diag_intstry(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Internal force term in force balance (y)
371
372      IF  ( iom_use( "normstr"  ) ) CALL iom_put( "normstr"     ,   diag_sig1(:,:)           * zswi(:,:) + zmiss(:,:) )   ! Normal stress
373      IF  ( iom_use( "sheastr"  ) ) CALL iom_put( "sheastr"     ,   diag_sig2(:,:)           * zswi(:,:) + zmiss(:,:) )   ! Shear stress
374
375      !--------------------------------
376      ! Global ice diagnostics (SIMIP)
377      !--------------------------------
378
379      IF ( iom_use( "NH_icearea" ) .OR. iom_use( "NH_icevolu" ) .OR. iom_use( "NH_iceextt" ) )   THEN   ! NH integrated diagnostics
380 
381         WHERE( ff_t > 0._wp ); zswi(:,:) = 1.0e-12
382         ELSEWHERE            ; zswi(:,:) = 0.
383         END WHERE
384
385         zdiag_area_nh = glob_sum( at_i(:,:) * zswi(:,:) * e1e2t(:,:) )
386         zdiag_volu_nh = glob_sum( vt_i(:,:) * zswi(:,:) * e1e2t(:,:) )
387
388         WHERE( ff_t > 0._wp .AND. at_i > 0.15 ); zswi(:,:) = 1.0e-12
389         ELSEWHERE                              ; zswi(:,:) = 0.
390         END WHERE
391
392         zdiag_extt_nh = glob_sum( zswi(:,:) * e1e2t(:,:) )
393
394         IF ( iom_use( "NH_icearea" ) ) CALL iom_put( "NH_icearea" ,  zdiag_area_nh  )
395         IF ( iom_use( "NH_icevolu" ) ) CALL iom_put( "NH_icevolu" ,  zdiag_volu_nh  )
396         IF ( iom_use( "NH_iceextt" ) ) CALL iom_put( "NH_iceextt" ,  zdiag_extt_nh  )
397
398      ENDIF
399
400      IF ( iom_use( "SH_icearea" ) .OR. iom_use( "SH_icevolu" ) .OR. iom_use( "SH_iceextt" ) )   THEN   ! SH integrated diagnostics
401
402         WHERE( ff_t < 0._wp ); zswi(:,:) = 1.0e-12; 
403         ELSEWHERE            ; zswi(:,:) = 0.
404         END WHERE
405
406         zdiag_area_sh = glob_sum( at_i(:,:) * zswi(:,:) * e1e2t(:,:) ) 
407         zdiag_volu_sh = glob_sum( vt_i(:,:) * zswi(:,:) * e1e2t(:,:) )
408
409         WHERE( ff_t < 0._wp .AND. at_i > 0.15 ); zswi(:,:) = 1.0e-12
410         ELSEWHERE                              ; zswi(:,:) = 0.
411         END WHERE
412
413         zdiag_extt_sh = glob_sum( zswi(:,:) * e1e2t(:,:) )
414
415         IF ( iom_use( "SH_icearea" ) ) CALL iom_put( "SH_icearea", zdiag_area_sh )
416         IF ( iom_use( "SH_icevolu" ) ) CALL iom_put( "SH_icevolu", zdiag_volu_sh )
417         IF ( iom_use( "SH_iceextt" ) ) CALL iom_put( "SH_iceextt", zdiag_extt_sh )
418
419      ENDIF 
420
421      !     !  Create an output files (output.lim.abort.nc) if S < 0 or u > 20 m/s
422      !     IF( kindic < 0 )   CALL ice_wri_state( 'output.abort' )
423      !     not yet implemented
424     
425      IF( nn_timing == 1 )  CALL timing_stop('icewri')
426     
427   END SUBROUTINE ice_wri
428
429 
430   SUBROUTINE ice_wri_state( kt, kid, kh_i )
431      !!---------------------------------------------------------------------
432      !!                 ***  ROUTINE ice_wri_state  ***
433      !!       
434      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
435      !!      the instantaneous ice state and forcing fields for ice model
436      !!        Used to find errors in the initial state or save the last
437      !!      ocean state in case of abnormal end of a simulation
438      !!
439      !! History :
440      !!   4.0  !  2013-06  (C. Rousset)
441      !!----------------------------------------------------------------------
442      INTEGER, INTENT( in )   ::   kt               ! ocean time-step index)
443      INTEGER, INTENT( in )   ::   kid , kh_i
444      INTEGER                 ::   nz_i, jl
445      REAL(wp), DIMENSION(jpl) :: jcat
446      !!----------------------------------------------------------------------
447      DO jl = 1, jpl
448         jcat(jl) = REAL(jl)
449      ENDDO
450     
451      CALL histvert( kid, "ncatice", "Ice Categories","", jpl, jcat, nz_i, "up")
452
453      CALL histdef( kid, "sithic", "Ice thickness"           , "m"      ,   &
454      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
455      CALL histdef( kid, "siconc", "Ice concentration"       , "%"      ,   &
456      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
457      CALL histdef( kid, "sitemp", "Ice temperature"         , "C"      ,   &
458      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
459      CALL histdef( kid, "sivelu", "i-Ice speed "            , "m/s"    ,   &
460      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
461      CALL histdef( kid, "sivelv", "j-Ice speed "            , "m/s"    ,   &
462      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
463      CALL histdef( kid, "sistru", "i-Wind stress over ice " , "Pa"     ,   &
464      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
465      CALL histdef( kid, "sistrv", "j-Wind stress over ice " , "Pa"     ,   &
466      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
467      CALL histdef( kid, "sisflx", "Solar flux over ocean"     , "w/m2" ,   &
468      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
469      CALL histdef( kid, "sinflx", "Non-solar flux over ocean" , "w/m2" ,   &
470      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
471      CALL histdef( kid, "isnowpre", "Snow precipitation"      , "kg/m2/s",   &
472      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
473      CALL histdef( kid, "sisali", "Ice salinity"            , "PSU"    ,   &
474      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
475      CALL histdef( kid, "sivolu", "Ice volume"              , "m"      ,   &
476      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
477      CALL histdef( kid, "sidive", "Ice divergence"          , "10-8s-1",   &
478      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
479
480      ! MV MP 2016
481      IF ( ln_pnd ) THEN
482         CALL histdef( kid, "si_amp", "Melt pond fraction"      , "%"      ,   &
483      &         jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
484         CALL histdef( kid, "si_vmp", "Melt pond volume"        ,  "m"     ,   &
485      &         jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
486      ENDIF
487      ! END MV MP 2016
488
489      CALL histdef( kid, "vfxbog", "Ice bottom production"   , "m/s"    ,   &
490      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
491      CALL histdef( kid, "vfxdyn", "Ice dynamic production"  , "m/s"    ,   &
492      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
493      CALL histdef( kid, "vfxopw", "Ice open water prod"     , "m/s"    ,   &
494      &       jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
495      CALL histdef( kid, "vfxsni", "Snow ice production "    , "m/s"    ,   &
496      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
497      CALL histdef( kid, "vfxres", "Ice prod from corrections" , "m/s"  ,   &
498      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
499      CALL histdef( kid, "vfxbom", "Ice bottom melt"         , "m/s"    ,   &
500      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
501      CALL histdef( kid, "vfxsum", "Ice surface melt"        , "m/s"    ,   &
502      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
503
504      CALL histdef( kid, "sithicat", "Ice thickness"         , "m"      ,   &
505      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
506      CALL histdef( kid, "siconcat", "Ice concentration"     , "%"      ,   &
507      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
508      CALL histdef( kid, "sisalcat", "Ice salinity"           , ""      ,   &
509      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
510      CALL histdef( kid, "sitemcat", "Ice temperature"       , "C"      ,   &
511      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
512      CALL histdef( kid, "snthicat", "Snw thickness"         , "m"      ,   &
513      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
514      CALL histdef( kid, "sntemcat", "Snw temperature"       , "C"      ,   &
515      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt )
516
517      CALL histend( kid, snc4set )   ! end of the file definition
518
519      CALL histwrite( kid, "sithic", kt, htm_i         , jpi*jpj, (/1/) )   
520      CALL histwrite( kid, "siconc", kt, at_i          , jpi*jpj, (/1/) )
521      CALL histwrite( kid, "sitemp", kt, tm_i - rt0    , jpi*jpj, (/1/) )
522      CALL histwrite( kid, "sivelu", kt, u_ice          , jpi*jpj, (/1/) )
523      CALL histwrite( kid, "sivelv", kt, v_ice          , jpi*jpj, (/1/) )
524      CALL histwrite( kid, "sistru", kt, utau_ice       , jpi*jpj, (/1/) )
525      CALL histwrite( kid, "sistrv", kt, vtau_ice       , jpi*jpj, (/1/) )
526      CALL histwrite( kid, "sisflx", kt, qsr , jpi*jpj, (/1/) )
527      CALL histwrite( kid, "sinflx", kt, qns , jpi*jpj, (/1/) )
528      CALL histwrite( kid, "isnowpre", kt, sprecip        , jpi*jpj, (/1/) )
529      CALL histwrite( kid, "sisali", kt, smt_i          , jpi*jpj, (/1/) )
530      CALL histwrite( kid, "sivolu", kt, vt_i           , jpi*jpj, (/1/) )
531      CALL histwrite( kid, "sidive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) )
532
533      ! MV MP 2016
534      IF ( ln_pnd ) THEN
535         CALL histwrite( kid, "si_amp", kt, at_ip         , jpi*jpj, (/1/) )
536         CALL histwrite( kid, "si_vmp", kt, vt_ip         , jpi*jpj, (/1/) )
537      ENDIF
538      ! END MV MP 2016
539
540      CALL histwrite( kid, "vfxbog", kt, wfx_bog        , jpi*jpj, (/1/) )
541      CALL histwrite( kid, "vfxdyn", kt, wfx_dyn        , jpi*jpj, (/1/) )
542      CALL histwrite( kid, "vfxopw", kt, wfx_opw        , jpi*jpj, (/1/) )
543      CALL histwrite( kid, "vfxsni", kt, wfx_sni        , jpi*jpj, (/1/) )
544      CALL histwrite( kid, "vfxres", kt, wfx_res        , jpi*jpj, (/1/) )
545      CALL histwrite( kid, "vfxbom", kt, wfx_bom        , jpi*jpj, (/1/) )
546      CALL histwrite( kid, "vfxsum", kt, wfx_sum        , jpi*jpj, (/1/) )
547      IF ( ln_pnd ) &
548         CALL histwrite( kid, "vfxpnd", kt, wfx_pnd     , jpi*jpj, (/1/) )
549
550      CALL histwrite( kid, "sithicat", kt, ht_i        , jpi*jpj*jpl, (/1/) )   
551      CALL histwrite( kid, "siconcat", kt, a_i         , jpi*jpj*jpl, (/1/) )   
552      CALL histwrite( kid, "sisalcat", kt, sm_i        , jpi*jpj*jpl, (/1/) )   
553      CALL histwrite( kid, "sitemcat", kt, tm_i - rt0  , jpi*jpj*jpl, (/1/) )   
554      CALL histwrite( kid, "snthicat", kt, ht_s        , jpi*jpj*jpl, (/1/) )   
555      CALL histwrite( kid, "sntemcat", kt, tm_su - rt0 , jpi*jpj*jpl, (/1/) )   
556
557      ! Close the file
558      ! -----------------
559      !CALL histclo( kid )
560
561    END SUBROUTINE ice_wri_state
562
563#else
564   !!----------------------------------------------------------------------
565   !!   Default option :         Empty module          NO LIM sea-ice model
566   !!----------------------------------------------------------------------
567CONTAINS
568   SUBROUTINE ice_wri          ! Empty routine
569   END SUBROUTINE ice_wri
570#endif
571
572   !!======================================================================
573END MODULE icewri
Note: See TracBrowser for help on using the repository browser.