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

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

compilation bugs

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