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.

icewri.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

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

Last change on this file since 8486 was 8486, checked in by clem, 7 years ago
changes in style - part1 - (now the code looks better txs to Gurvan's comments)
File size: 36.9 KB

Rev	Line
[8413]	1	MODULE 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	!!----------------------------------------------------------------------
[8486]	10	!! ice_wri : write of the diagnostics variables in ouput file
[8413]	11	!! ice_wri_state : write for initial state or/and abandon
	12	!!----------------------------------------------------------------------
[8486]	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
[8413]	20	!
[8486]	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
[8413]	28
	29	IMPLICIT NONE
	30	PRIVATE
	31
[8486]	32	PUBLIC ice_wri ! called by ice_stp
[8413]	33	PUBLIC ice_wri_state ! called by dia_wri_state
	34
	35	!!----------------------------------------------------------------------
[8486]	36	!! NEMO/ICE 4.0 , NEMO Consortium (2017)
[8413]	37	!! $Id: icewri.F90 8409 2017-08-07 15:29:21Z clem $
	38	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	39	!!----------------------------------------------------------------------
	40	CONTAINS
	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
[8486]	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
[8413]	61	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zswi2, zmiss2
[8486]	62	!
[8413]	63	! Global ice diagnostics (SIMIP)
[8486]	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
[8413]	66	!!-------------------------------------------------------------------
	67
[8486]	68	IF( nn_timing == 1 ) CALL timing_start('icewri')
[8413]	69
	70	!----------------------------------------
	71	! Brine volume, switches, missing values
	72	!----------------------------------------
	73
[8486]	74	CALL ice_var_bv ! brine volume
[8413]	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
[8486]	92	zmiss_val = 1.e20_wp
	93	zmiss (:,:) = zmiss_val * ( 1._wp - zswi (:,:) )
	94	zmiss2(:,:,:) = zmiss_val * ( 1._wp - zswi2(:,:,:) )
[8413]	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(:,:) )
[8486]	108	!!gm I don't understand the variable below.... why not multiplied by a_i_b or (1-a_i_b) ???
[8413]	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(:,:) )
[8486]	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)
[8413]	113
	114	! velocity
[8486]	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
[8413]	117
[8486]	118	IF( iom_use('icevel') .OR. iom_use('icevel_mv') ) THEN
[8413]	119	DO jj = 2 , jpjm1
	120	DO ji = 2 , jpim1
[8486]	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) = O.5_wp * SQRT( z2da * z2da + z2db * z2db )
[8413]	124	END DO
	125	END DO
	126	CALL lbc_lnk( z2d, 'T', 1. )
[8486]	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)
[8413]	129	ENDIF
	130
[8486]	131	IF( iom_use('tau_icebfr') ) CALL iom_put( "tau_icebfr" , tau_icebfr ) ! ice friction with ocean bottom (landfast ice)
[8413]	132	!
[8486]	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
[8413]	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
[8486]	146	CALL iom_put( "snowpre" , sprecip * 86400. ) ! snow precipitation [m/day]
[8413]	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
[8486]	221	!!gm ====>>>>> THIS should be moved in icerhg_evp (generalize this everywhere it is possible and logic...)
[8413]	222	! specific outputs for EVP rheology
[8486]	223	IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN
[8413]	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
[8486]	249	!!gm <<<<<<======= end
[8413]	250
[8486]	251	!!gm ====>>>>> THIS should be moved where at_ip, vt_ip are computed fro the last time in the time-step (limmpd)
[8413]	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
[8486]	258	!!gm <<<<<<======= end
[8413]	259
	260	!----------------------------------
	261	! Output category-dependent fields
	262	!----------------------------------
[8486]	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
[8413]	271
	272	! MV MP 2016
	273	IF ( ln_pnd ) THEN
[8486]	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
[8413]	278	ENDIF
	279	! END MV MP 2016
	280
	281	!--------------------------------
	282	! Add-ons for SIMIP
	283	!--------------------------------
	284	zrho1 = ( rau0 - rhoic ) / rau0; zrho2 = rhosn / rau0
	285
[8486]	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
[8413]	291
[8486]	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
[8413]	295
[8486]	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
[8413]	302
[8486]	303	IF ( iom_use( "icefb" ) ) THEN
[8413]	304	zfb(:,:) = ( zrho1 * htm_i(:,:) - zrho2 * htm_s(:,:) )
	305	WHERE( zfb < 0._wp ) ; zfb = 0._wp ; END WHERE
[8486]	306	CALL iom_put( "icefb" , zfb(:,:) * zswi(:,:) + zmiss(:,:) ) ! Ice freeboard
[8413]	307	ENDIF
	308
[8486]	309	IF ( iom_use( "isnhcneg" ) ) CALL iom_put( "isnhcneg" , - et_s(:,:) * zswi(:,:) + zmiss(:,:) ) ! Snow total heat content
[8413]	310
[8486]	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
[8413]	319
[8486]	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
[8413]	323
[8486]	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
[8413]	326
[8486]	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
[8413]	329
[8486]	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
[8413]	333
[8486]	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
[8413]	337
[8486]	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
[8413]	340
[8486]	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
[8413]	343
[8486]	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
[8413]	346
[8486]	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
[8413]	349
[8486]	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)
[8413]	352
[8486]	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)
[8413]	355
[8486]	356	IF ( iom_use( "icestr" ) ) CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) + zmiss(:,:) ) ! Ice strength
[8413]	357
[8486]	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)
[8413]	360
[8486]	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)
[8413]	363
[8486]	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)
[8413]	366
[8486]	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
[8413]	369
	370	!--------------------------------
	371	! Global ice diagnostics (SIMIP)
	372	!--------------------------------
[8486]	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.
[8413]	378	END WHERE
	379	zdiag_area_nh = glob_sum( at_i(:,:) * zswi(:,:) * e1e2t(:,:) )
	380	zdiag_volu_nh = glob_sum( vt_i(:,:) * zswi(:,:) * e1e2t(:,:) )
[8486]	381	!
	382	WHERE( ff_t > 0._wp .AND. at_i > 0.15 ) ; zswi(:,:) = 1.0e-12
	383	ELSEWHERE ; zswi(:,:) = 0.
[8413]	384	END WHERE
	385	zdiag_extt_nh = glob_sum( zswi(:,:) * e1e2t(:,:) )
[8486]	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	!
[8413]	391	ENDIF
[8486]	392	!
	393	IF ( iom_use( "SH_icearea" ) .OR. iom_use( "SH_icevolu" ) .OR. iom_use( "SH_iceextt" ) ) THEN ! SH diagnostics
	394	!
[8413]	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(:,:) )
[8486]	400	!
[8413]	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(:,:) )
[8486]	405	!
[8413]	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 )
[8486]	409	!
[8413]	410	ENDIF
[8486]	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	!
[8413]	417	IF( nn_timing == 1 ) CALL timing_stop('icewri')
[8486]	418	!
[8413]	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	!!
[8486]	431	!! History : 4.0 ! 2013-06 (C. Rousset)
[8413]	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	!!----------------------------------------------------------------------
[8486]	438	!
[8413]	439	DO jl = 1, jpl
	440	jcat(jl) = REAL(jl)
[8486]	441	END DO
[8413]	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_i1.0e8 , jpijpj, (/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 , jpijpjjpl, (/1/) )
	543	CALL histwrite( kid, "siconcat", kt, a_i , jpijpjjpl, (/1/) )
	544	CALL histwrite( kid, "sisalcat", kt, sm_i , jpijpjjpl, (/1/) )
	545	CALL histwrite( kid, "sitemcat", kt, tm_i - rt0 , jpijpjjpl, (/1/) )
	546	CALL histwrite( kid, "snthicat", kt, ht_s , jpijpjjpl, (/1/) )
	547	CALL histwrite( kid, "sntemcat", kt, tm_su - rt0 , jpijpjjpl, (/1/) )
	548
	549	! Close the file
	550	! -----------------
[8486]	551	!!gm I don't understand why the file is not closed !
[8413]	552	!CALL histclo( kid )
[8486]	553	!
[8413]	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	!!======================================================================
	563	END MODULE icewri

Note: See TracBrowser for help on using the repository browser.

Download in other formats: