New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
limwri.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90 @ 8316

Last change on this file since 8316 was 8316, checked in by clem, 7 years ago

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