Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

iceistate.F90 @ 12636

Last change on this file since 12636 was 12636, checked in by dancopsey, 4 years ago

Add:

Melt pond lids
Melt pond maximum area and thickness
Melt pond vertical flushing
Area contributing to melt ponds depends on total ice fraction

File size: 27.1 KB

Rev	Line
[8586]	1	MODULE iceistate
	2	!!======================================================================
	3	!! * MODULE iceistate *
	4	!! sea-ice : Initialization of ice variables
	5	!!======================================================================
[9604]	6	!! History : 2.0 ! 2004-01 (C. Ethe, G. Madec) Original code
	7	!! 3.0 ! 2007 (M. Vancoppenolle) Rewrite for ice cats
	8	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
[8586]	9	!!----------------------------------------------------------------------
[9570]	10	#if defined key_si3
[8586]	11	!!----------------------------------------------------------------------
[9570]	12	!! 'key_si3' SI3 sea-ice model
[8586]	13	!!----------------------------------------------------------------------
	14	!! ice_istate : initialization of diagnostics ice variables
	15	!! ice_istate_init : initialization of ice state and namelist read
	16	!!----------------------------------------------------------------------
	17	USE phycst ! physical constant
	18	USE oce ! dynamics and tracers variables
	19	USE dom_oce ! ocean domain
	20	USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd
	21	USE sbc_ice , ONLY : tn_ice, snwice_mass, snwice_mass_b
	22	USE eosbn2 ! equation of state
	23	USE domvvl ! Variable volume
[11536]	24	USE ice ! sea-ice: variables
	25	USE ice1D ! sea-ice: thermodynamics variables
	26	USE icetab ! sea-ice: 1D <==> 2D transformation
	27	USE icevar ! sea-ice: operations
[8586]	28	!
	29	USE in_out_manager ! I/O manager
	30	USE iom ! I/O manager library
	31	USE lib_mpp ! MPP library
	32	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
	33	USE fldread ! read input fields
	34
	35	IMPLICIT NONE
	36	PRIVATE
	37
	38	PUBLIC ice_istate ! called by icestp.F90
	39	PUBLIC ice_istate_init ! called by icestp.F90
	40	!
[8813]	41	! !! namelist (namini)
[11536]	42	LOGICAL, PUBLIC :: ln_iceini !: Ice initialization or not
	43	LOGICAL, PUBLIC :: ln_iceini_file !: Ice initialization from 2D netcdf file
	44	REAL(wp) :: rn_thres_sst
	45	REAL(wp) :: rn_hti_ini_n, rn_hts_ini_n, rn_ati_ini_n, rn_smi_ini_n, rn_tmi_ini_n, rn_tsu_ini_n, rn_tms_ini_n
	46	REAL(wp) :: rn_hti_ini_s, rn_hts_ini_s, rn_ati_ini_s, rn_smi_ini_s, rn_tmi_ini_s, rn_tsu_ini_s, rn_tms_ini_s
	47	REAL(wp) :: rn_apd_ini_n, rn_hpd_ini_n
	48	REAL(wp) :: rn_apd_ini_s, rn_hpd_ini_s
	49	!
	50	! ! if ln_iceini_file = T
	51	INTEGER , PARAMETER :: jpfldi = 9 ! maximum number of files to read
	52	INTEGER , PARAMETER :: jp_hti = 1 ! index of ice thickness (m)
	53	INTEGER , PARAMETER :: jp_hts = 2 ! index of snw thickness (m)
	54	INTEGER , PARAMETER :: jp_ati = 3 ! index of ice fraction (-)
	55	INTEGER , PARAMETER :: jp_smi = 4 ! index of ice salinity (g/kg)
	56	INTEGER , PARAMETER :: jp_tmi = 5 ! index of ice temperature (K)
	57	INTEGER , PARAMETER :: jp_tsu = 6 ! index of ice surface temp (K)
	58	INTEGER , PARAMETER :: jp_tms = 7 ! index of snw temperature (K)
	59	INTEGER , PARAMETER :: jp_apd = 8 ! index of pnd fraction (-)
	60	INTEGER , PARAMETER :: jp_hpd = 9 ! index of pnd depth (m)
	61	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: si ! structure of input fields (file informations, fields read)
	62	!
[8586]	63	!!----------------------------------------------------------------------
[9598]	64	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	65	!! $Id$
[10413]	66	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[8586]	67	!!----------------------------------------------------------------------
	68	CONTAINS
	69
[11536]	70	SUBROUTINE ice_istate( kt )
[8586]	71	!!-------------------------------------------------------------------
	72	!! * ROUTINE ice_istate *
	73	!!
	74	!! ** Purpose : defined the sea-ice initial state
	75	!!
	76	!! ** Method : This routine will put some ice where ocean
	77	!! is at the freezing point, then fill in ice
	78	!! state variables using prescribed initial
	79	!! values in the namelist
	80	!!
	81	!! ** Steps : 1) Set initial surface and basal temperatures
	82	!! 2) Recompute or read sea ice state variables
	83	!! 3) Fill in the ice thickness distribution using gaussian
	84	!! 4) Fill in space-dependent arrays for state variables
	85	!! 5) snow-ice mass computation
	86	!! 6) store before fields
	87	!!
	88	!! ** Notes : o_i, t_su, t_s, t_i, sz_i must be filled everywhere, even
[11536]	89	!! where there is no ice
[8586]	90	!!--------------------------------------------------------------------
[11536]	91	INTEGER, INTENT(in) :: kt ! time step
	92	!!
[8813]	93	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[11536]	94	REAL(wp) :: ztmelts
[8586]	95	INTEGER , DIMENSION(4) :: itest
	96	REAL(wp), DIMENSION(jpi,jpj) :: z2d
	97	REAL(wp), DIMENSION(jpi,jpj) :: zswitch ! ice indicator
[11536]	98	REAL(wp), DIMENSION(jpi,jpj) :: zht_i_ini, zat_i_ini, ztm_s_ini !data from namelist or nc file
	99	REAL(wp), DIMENSION(jpi,jpj) :: zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
	100	REAL(wp), DIMENSION(jpi,jpj) :: zapnd_ini, zhpnd_ini !data from namelist or nc file
	101	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zti_3d , zts_3d !temporary arrays
	102	!!
	103	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d
[8586]	104	!--------------------------------------------------------------------
	105
	106	IF(lwp) WRITE(numout,*)
	107	IF(lwp) WRITE(numout,*) 'ice_istate: sea-ice initialization '
	108	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
	109
[11536]	110	!---------------------------
	111	! 1) 1st init. of the fields
	112	!---------------------------
[8586]	113	!
[11536]	114	! basal temperature (considered at freezing point) [Kelvin]
	115	CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
	116	t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)
	117	!
	118	! surface temperature and conductivity
[8586]	119	DO jl = 1, jpl
[9416]	120	t_su (:,:,jl) = rt0 * tmask(:,:,1) ! temp at the surface
[10534]	121	cnd_ice(:,:,jl) = 0._wp ! initialisation of the effective conductivity at the top of ice/snow (ln_cndflx=T)
[8586]	122	END DO
[8813]	123	!
[11536]	124	! ice and snw temperatures
	125	DO jl = 1, jpl
	126	DO jk = 1, nlay_i
	127	t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
	128	END DO
	129	DO jk = 1, nlay_s
	130	t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
	131	END DO
	132	END DO
	133	!
	134	! specific temperatures for coupled runs
	135	tn_ice (:,:,:) = t_i (:,:,1,:)
	136	t1_ice (:,:,:) = t_i (:,:,1,:)
[8586]	137
[11536]	138	! heat contents
	139	e_i (:,:,:,:) = 0._wp
	140	e_s (:,:,:,:) = 0._wp
	141
	142	! general fields
	143	a_i (:,:,:) = 0._wp
	144	v_i (:,:,:) = 0._wp
	145	v_s (:,:,:) = 0._wp
	146	sv_i(:,:,:) = 0._wp
	147	oa_i(:,:,:) = 0._wp
	148	!
	149	h_i (:,:,:) = 0._wp
	150	h_s (:,:,:) = 0._wp
	151	s_i (:,:,:) = 0._wp
	152	o_i (:,:,:) = 0._wp
	153	!
	154	! melt ponds
	155	a_ip (:,:,:) = 0._wp
	156	v_ip (:,:,:) = 0._wp
[12636]	157	lh_ip (:,:,:) = 0._wp
[11536]	158	a_ip_frac(:,:,:) = 0._wp
[12636]	159	a_ip_eff (:,:,:) = 0._wp
[11536]	160	h_ip (:,:,:) = 0._wp
	161	!
	162	! ice velocities
	163	u_ice (:,:) = 0._wp
	164	v_ice (:,:) = 0._wp
	165	!
	166	!------------------------------------------------------------------------
	167	! 2) overwrite some of the fields with namelist parameters or netcdf file
	168	!------------------------------------------------------------------------
[8586]	169	IF( ln_iceini ) THEN
	170	! !---------------!
	171	IF( ln_iceini_file )THEN ! Read a file !
	172	! !---------------!
[11536]	173	WHERE( ff_t(:,:) >= 0._wp ) ; zswitch(:,:) = 1._wp
	174	ELSEWHERE ; zswitch(:,:) = 0._wp
	175	END WHERE
[8586]	176	!
[11536]	177	CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step
[8586]	178	!
[11536]	179	! -- mandatory fields -- !
	180	zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1)
	181	zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1)
	182	zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1)
	183
	184	! -- optional fields -- !
	185	! if fields do not exist then set them to the values present in the namelist (except for snow and surface temperature)
	186	!
	187	! ice salinity
	188	IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) &
	189	& si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
	190	zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1)
	191	!
	192	! ice temperature
	193	IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) &
	194	& si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
	195	ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1)
	196	!
	197	! surface temperature => set to ice temperature if it exists
	198	IF ( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN
	199	si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
	200	ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN
	201	si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
	202	ENDIF
	203	zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1)
	204	!
	205	! snow temperature => set to ice temperature if it exists
	206	IF ( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN
	207	si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
	208	ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN
	209	si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
	210	ENDIF
	211	ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1)
	212	!
	213	! pond concentration
	214	IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) &
	215	& si(jp_apd)%fnow(:,:,1) = ( rn_apd_ini_n * zswitch + rn_apd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) & ! rn_apd = pond fraction => rn_apnd * a_i = pond conc.
	216	& * si(jp_ati)%fnow(:,:,1)
	217	zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1)
	218	!
	219	! pond depth
	220	IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) &
	221	& si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
	222	zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1)
	223	!
	224	! change the switch for the following
	225	WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1)
	226	ELSEWHERE ; zswitch(:,:) = 0._wp
[8586]	227	END WHERE
[8813]	228	! !---------------!
[8586]	229	ELSE ! Read namelist !
	230	! !---------------!
[11536]	231	! no ice if (sst - Tfreez) >= thresold
[8813]	232	WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp
	233	ELSEWHERE ; zswitch(:,:) = tmask(:,:,1)
[8586]	234	END WHERE
[8813]	235	!
[8586]	236	! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
	237	WHERE( ff_t(:,:) >= 0._wp )
	238	zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
	239	zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
	240	zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
	241	zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
	242	ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
[11536]	243	zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:)
	244	ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:)
	245	zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
	246	zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:)
[8586]	247	ELSEWHERE
	248	zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
	249	zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
	250	zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
	251	zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
	252	ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
[11536]	253	zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:)
	254	ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:)
	255	zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
	256	zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:)
[8586]	257	END WHERE
	258	!
	259	ENDIF
[11536]	260
	261	! make sure ponds = 0 if no ponds scheme
	262	IF ( .NOT.ln_pnd ) THEN
	263	zapnd_ini(:,:) = 0._wp
	264	zhpnd_ini(:,:) = 0._wp
	265	ENDIF
[8586]	266
[11536]	267	!-------------!
	268	! fill fields !
	269	!-------------!
	270	! select ice covered grid points
	271	npti = 0 ; nptidx(:) = 0
	272	DO jj = 1, jpj
	273	DO ji = 1, jpi
	274	IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
	275	npti = npti + 1
	276	nptidx(npti) = (jj - 1) * jpi + ji
	277	ENDIF
	278	END DO
	279	END DO
[8586]	280
[11536]	281	! move to 1D arrays: (jpi,jpj) -> (jpi*jpj)
	282	CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti) , zht_i_ini )
	283	CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti) , zht_s_ini )
	284	CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti) , zat_i_ini )
	285	CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini )
	286	CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini )
	287	CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti) , zt_su_ini )
	288	CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti) , zsm_i_ini )
	289	CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti) , zapnd_ini )
	290	CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti) , zhpnd_ini )
[10413]	291
[11536]	292	! allocate temporary arrays
	293	ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), &
	294	& zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) )
	295
	296	! distribute 1-cat into jpl-cat: (jpijpj) -> (jpijpj,jpl)
	297	CALL ice_var_itd( h_i_1d(1:npti) , h_s_1d(1:npti) , at_i_1d(1:npti), &
	298	& zhi_2d , zhs_2d , zai_2d , &
	299	& t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti), a_ip_1d(1:npti), h_ip_1d(1:npti), &
	300	& zti_2d , zts_2d , ztsu_2d , zsi_2d , zaip_2d , zhip_2d )
[8586]	301
[11536]	302	! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
	303	DO jl = 1, jpl
	304	zti_3d(:,:,jl) = rt0 * tmask(:,:,1)
	305	zts_3d(:,:,jl) = rt0 * tmask(:,:,1)
	306	END DO
	307	CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d , h_i )
	308	CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d , h_s )
	309	CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d , a_i )
	310	CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d , zti_3d )
	311	CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d , zts_3d )
	312	CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d , t_su )
	313	CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d , s_i )
	314	CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d , a_ip )
	315	CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d , h_ip )
[8586]	316
[11536]	317	! deallocate temporary arrays
	318	DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
	319	& zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d )
	320
	321	! calculate extensive and intensive variables
	322	CALL ice_var_salprof ! for sz_i
	323	DO jl = 1, jpl
[8586]	324	DO jj = 1, jpj
	325	DO ji = 1, jpi
[11536]	326	v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
	327	v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
	328	sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl)
[8586]	329	END DO
	330	END DO
	331	END DO
[8813]	332	!
[11536]	333	DO jl = 1, jpl
	334	DO jk = 1, nlay_s
[8586]	335	DO jj = 1, jpj
	336	DO ji = 1, jpi
[11536]	337	t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
	338	e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
	339	& rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
[8586]	340	END DO
	341	END DO
	342	END DO
	343	END DO
[8813]	344	!
[11536]	345	DO jl = 1, jpl
	346	DO jk = 1, nlay_i
[8586]	347	DO jj = 1, jpj
	348	DO ji = 1, jpi
[11536]	349	t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl)
	350	ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
	351	e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
	352	& rhoi * ( rcpi * ( ztmelts - t_i(ji,jj,jk,jl) ) + &
	353	& rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) &
	354	& - rcp * ( ztmelts - rt0 ) )
[8586]	355	END DO
	356	END DO
	357	END DO
	358	END DO
[8933]	359
[11536]	360	! Melt ponds
	361	WHERE( a_i > epsi10 )
	362	a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
	363	ELSEWHERE
	364	a_ip_frac(:,:,:) = 0._wp
	365	END WHERE
[12636]	366	a_ip_eff(:,:,:) = a_ip_frac(:,:,:)
[11536]	367	v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
	368
	369	! specific temperatures for coupled runs
	370	tn_ice(:,:,:) = t_su(:,:,:)
	371	t1_ice(:,:,:) = t_i (:,:,1,:)
[8813]	372	!
[8586]	373	ENDIF ! ln_iceini
[8813]	374	!
[11536]	375	at_i(:,:) = SUM( a_i, dim=3 )
[8586]	376	!
	377	!----------------------------------------------
[11536]	378	! 3) Snow-ice mass (case ice is fully embedded)
[8586]	379	!----------------------------------------------
[9935]	380	snwice_mass (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3 ) ! snow+ice mass
[8586]	381	snwice_mass_b(:,:) = snwice_mass(:,:)
	382	!
	383	IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area
[8813]	384	!
[8586]	385	sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
	386	sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
[8813]	387	!
[8586]	388	IF( .NOT.ln_linssh ) THEN
[8813]	389	!
[8586]	390	WHERE( ht_0(:,:) > 0 ) ; z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:)
	391	ELSEWHERE ; z2d(:,:) = 1._wp ; END WHERE
[8813]	392	!
[8586]	393	DO jk = 1,jpkm1 ! adjust initial vertical scale factors
	394	e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:)
	395	e3t_b(:,:,jk) = e3t_n(:,:,jk)
	396	e3t_a(:,:,jk) = e3t_n(:,:,jk)
	397	END DO
[8813]	398	!
[8586]	399	! Reconstruction of all vertical scale factors at now and before time-steps
	400	! =========================================================================
	401	! Horizontal scale factor interpolations
	402	! --------------------------------------
	403	CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
	404	CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
	405	CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
	406	CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
	407	CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
	408	! Vertical scale factor interpolations
	409	! ------------------------------------
	410	CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' )
	411	CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
	412	CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
	413	CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
	414	CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
	415	! t- and w- points depth
	416	! ----------------------
	417	!!gm not sure of that....
	418	gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
	419	gdepw_n(:,:,1) = 0.0_wp
	420	gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
	421	DO jk = 2, jpk
	422	gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk )
	423	gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
	424	gde3w_n(:,:,jk) = gdept_n(:,:,jk ) - sshn (:,:)
	425	END DO
	426	ENDIF
	427	ENDIF
	428
	429	!------------------------------------
[11536]	430	! 4) store fields at before time-step
[8586]	431	!------------------------------------
	432	! it is only necessary for the 1st interpolation by Agrif
	433	a_i_b (:,:,:) = a_i (:,:,:)
	434	e_i_b (:,:,:,:) = e_i (:,:,:,:)
	435	v_i_b (:,:,:) = v_i (:,:,:)
	436	v_s_b (:,:,:) = v_s (:,:,:)
	437	e_s_b (:,:,:,:) = e_s (:,:,:,:)
	438	sv_i_b (:,:,:) = sv_i (:,:,:)
	439	oa_i_b (:,:,:) = oa_i (:,:,:)
	440	u_ice_b(:,:) = u_ice(:,:)
	441	v_ice_b(:,:) = v_ice(:,:)
[10527]	442	! total concentration is needed for Lupkes parameterizations
	443	at_i_b (:,:) = at_i (:,:)
[8586]	444
[8885]	445	!!clem: output of initial state should be written here but it is impossible because
	446	!! the ocean and ice are in the same file
[10425]	447	!! CALL dia_wri_state( 'output.init' )
[8813]	448	!
[8586]	449	END SUBROUTINE ice_istate
	450
[8813]	451
[8586]	452	SUBROUTINE ice_istate_init
	453	!!-------------------------------------------------------------------
	454	!! * ROUTINE ice_istate_init *
	455	!!
[8813]	456	!! ** Purpose : Definition of initial state of the ice
[8586]	457	!!
[8813]	458	!! ** Method : Read the namini namelist and check the parameter
	459	!! values called at the first timestep (nit000)
[8586]	460	!!
[8813]	461	!! ** input : Namelist namini
[8586]	462	!!
	463	!!-----------------------------------------------------------------------------
[11229]	464	INTEGER :: ios ! Local integer output status for namelist read
[8813]	465	INTEGER :: ifpr, ierror
[8586]	466	!
	467	CHARACTER(len=256) :: cn_dir ! Root directory for location of ice files
[11536]	468	TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd
[8586]	469	TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read
	470	!
[11536]	471	NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, &
	472	& rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, &
	473	& rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, &
	474	& rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, &
	475	& rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, &
	476	& sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, cn_dir
[8586]	477	!!-----------------------------------------------------------------------------
	478	!
	479	REWIND( numnam_ice_ref ) ! Namelist namini in reference namelist : Ice initial state
	480	READ ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
[11536]	481	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist' )
[8586]	482	REWIND( numnam_ice_cfg ) ! Namelist namini in configuration namelist : Ice initial state
	483	READ ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
[11536]	484	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist' )
[8586]	485	IF(lwm) WRITE ( numoni, namini )
[8813]	486	!
[8586]	487	slf_i(jp_hti) = sn_hti ; slf_i(jp_hts) = sn_hts
[11536]	488	slf_i(jp_ati) = sn_ati ; slf_i(jp_smi) = sn_smi
	489	slf_i(jp_tmi) = sn_tmi ; slf_i(jp_tsu) = sn_tsu ; slf_i(jp_tms) = sn_tms
	490	slf_i(jp_apd) = sn_apd ; slf_i(jp_hpd) = sn_hpd
[8586]	491	!
	492	IF(lwp) THEN ! control print
	493	WRITE(numout,*)
	494	WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
	495	WRITE(numout,*) '~~~~~~~~~~~~~~~'
	496	WRITE(numout,*) ' Namelist namini:'
[11536]	497	WRITE(numout,*) ' ice initialization (T) or not (F) ln_iceini = ', ln_iceini
	498	WRITE(numout,*) ' ice initialization from a netcdf file ln_iceini_file = ', ln_iceini_file
	499	WRITE(numout,*) ' max ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst
	500	IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN
	501	WRITE(numout,*) ' initial snw thickness in the north-south rn_hts_ini = ', rn_hts_ini_n,rn_hts_ini_s
	502	WRITE(numout,*) ' initial ice thickness in the north-south rn_hti_ini = ', rn_hti_ini_n,rn_hti_ini_s
	503	WRITE(numout,*) ' initial ice concentr in the north-south rn_ati_ini = ', rn_ati_ini_n,rn_ati_ini_s
	504	WRITE(numout,*) ' initial ice salinity in the north-south rn_smi_ini = ', rn_smi_ini_n,rn_smi_ini_s
	505	WRITE(numout,*) ' initial surf temperat in the north-south rn_tsu_ini = ', rn_tsu_ini_n,rn_tsu_ini_s
	506	WRITE(numout,*) ' initial ice temperat in the north-south rn_tmi_ini = ', rn_tmi_ini_n,rn_tmi_ini_s
	507	WRITE(numout,*) ' initial snw temperat in the north-south rn_tms_ini = ', rn_tms_ini_n,rn_tms_ini_s
	508	WRITE(numout,*) ' initial pnd fraction in the north-south rn_apd_ini = ', rn_apd_ini_n,rn_apd_ini_s
	509	WRITE(numout,*) ' initial pnd depth in the north-south rn_hpd_ini = ', rn_hpd_ini_n,rn_hpd_ini_s
	510	ENDIF
[8586]	511	ENDIF
[8813]	512	!
[8586]	513	IF( ln_iceini_file ) THEN ! Ice initialization using input file
	514	!
	515	! set si structure
	516	ALLOCATE( si(jpfldi), STAT=ierror )
	517	IF( ierror > 0 ) THEN
[11536]	518	CALL ctl_stop( 'ice_istate_ini in iceistate: unable to allocate si structure' ) ; RETURN
[8586]	519	ENDIF
[8813]	520	!
[8586]	521	DO ifpr = 1, jpfldi
	522	ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
[11536]	523	IF( slf_i(ifpr)%ln_tint ) ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
[8586]	524	END DO
[8813]	525	!
[8586]	526	! fill si with slf_i and control print
[11536]	527	CALL fld_fill( si, slf_i, cn_dir, 'ice_istate_ini', 'initialization of sea ice fields', 'numnam_ice' )
[8813]	528	!
[8586]	529	ENDIF
[8813]	530	!
[11536]	531	IF( .NOT.ln_pnd ) THEN
	532	rn_apd_ini_n = 0. ; rn_apd_ini_s = 0.
	533	rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0.
	534	CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 when no ponds' )
	535	ENDIF
	536	!
[8586]	537	END SUBROUTINE ice_istate_init
	538
	539	#else
	540	!!----------------------------------------------------------------------
[9570]	541	!! Default option : Empty module NO SI3 sea-ice model
[8586]	542	!!----------------------------------------------------------------------
	543	#endif
	544
	545	!!======================================================================
	546	END MODULE iceistate

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceistate.F90 @ 12636

Download in other formats: