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.

iceistate.F90 in NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/iceistate.F90 @ 11927

Last change on this file since 11927 was 11927, checked in by dancopsey, 4 years ago
Fix compile errors and resolve conflicts
File size: 27.3 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
[11918]	43	INTEGER, PUBLIC :: nn_iceini_file ! Ice initialization:
	44	! 0 = Initialise sea ice based on SSTs
	45	! 1 = Initialise sea ice from single category netcdf file
	46	! 2 = Initialise sea ice from multi category restart file
[11536]	47	REAL(wp) :: rn_thres_sst
	48	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
	49	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
	50	REAL(wp) :: rn_apd_ini_n, rn_hpd_ini_n
	51	REAL(wp) :: rn_apd_ini_s, rn_hpd_ini_s
	52	!
[11927]	53	! ! if nn_iceini_file = 1
[11536]	54	INTEGER , PARAMETER :: jpfldi = 9 ! maximum number of files to read
	55	INTEGER , PARAMETER :: jp_hti = 1 ! index of ice thickness (m)
	56	INTEGER , PARAMETER :: jp_hts = 2 ! index of snw thickness (m)
	57	INTEGER , PARAMETER :: jp_ati = 3 ! index of ice fraction (-)
	58	INTEGER , PARAMETER :: jp_smi = 4 ! index of ice salinity (g/kg)
	59	INTEGER , PARAMETER :: jp_tmi = 5 ! index of ice temperature (K)
	60	INTEGER , PARAMETER :: jp_tsu = 6 ! index of ice surface temp (K)
	61	INTEGER , PARAMETER :: jp_tms = 7 ! index of snw temperature (K)
	62	INTEGER , PARAMETER :: jp_apd = 8 ! index of pnd fraction (-)
	63	INTEGER , PARAMETER :: jp_hpd = 9 ! index of pnd depth (m)
	64	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: si ! structure of input fields (file informations, fields read)
	65	!
[8586]	66	!!----------------------------------------------------------------------
[9598]	67	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	68	!! $Id$
[10413]	69	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[8586]	70	!!----------------------------------------------------------------------
	71	CONTAINS
	72
[11536]	73	SUBROUTINE ice_istate( kt )
[8586]	74	!!-------------------------------------------------------------------
	75	!! * ROUTINE ice_istate *
	76	!!
	77	!! ** Purpose : defined the sea-ice initial state
	78	!!
	79	!! ** Method : This routine will put some ice where ocean
	80	!! is at the freezing point, then fill in ice
	81	!! state variables using prescribed initial
	82	!! values in the namelist
	83	!!
	84	!! ** Steps : 1) Set initial surface and basal temperatures
	85	!! 2) Recompute or read sea ice state variables
	86	!! 3) Fill in the ice thickness distribution using gaussian
	87	!! 4) Fill in space-dependent arrays for state variables
	88	!! 5) snow-ice mass computation
	89	!! 6) store before fields
	90	!!
	91	!! ** Notes : o_i, t_su, t_s, t_i, sz_i must be filled everywhere, even
[11536]	92	!! where there is no ice
[8586]	93	!!--------------------------------------------------------------------
[11536]	94	INTEGER, INTENT(in) :: kt ! time step
	95	!!
[8813]	96	INTEGER :: ji, jj, jk, jl ! dummy loop indices
[11536]	97	REAL(wp) :: ztmelts
[8586]	98	INTEGER , DIMENSION(4) :: itest
	99	REAL(wp), DIMENSION(jpi,jpj) :: z2d
	100	REAL(wp), DIMENSION(jpi,jpj) :: zswitch ! ice indicator
[11536]	101	REAL(wp), DIMENSION(jpi,jpj) :: zht_i_ini, zat_i_ini, ztm_s_ini !data from namelist or nc file
	102	REAL(wp), DIMENSION(jpi,jpj) :: zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
	103	REAL(wp), DIMENSION(jpi,jpj) :: zapnd_ini, zhpnd_ini !data from namelist or nc file
	104	REAL(wp), DIMENSION(jpi,jpj,jpl) :: zti_3d , zts_3d !temporary arrays
	105	!!
	106	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d
[8586]	107	!--------------------------------------------------------------------
	108
	109	IF(lwp) WRITE(numout,*)
	110	IF(lwp) WRITE(numout,*) 'ice_istate: sea-ice initialization '
	111	IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
	112
[11536]	113	!---------------------------
	114	! 1) 1st init. of the fields
	115	!---------------------------
[8586]	116	!
[11536]	117	! basal temperature (considered at freezing point) [Kelvin]
	118	CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
	119	t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)
	120	!
	121	! surface temperature and conductivity
[8586]	122	DO jl = 1, jpl
[9416]	123	t_su (:,:,jl) = rt0 * tmask(:,:,1) ! temp at the surface
[10534]	124	cnd_ice(:,:,jl) = 0._wp ! initialisation of the effective conductivity at the top of ice/snow (ln_cndflx=T)
[8586]	125	END DO
[8813]	126	!
[11536]	127	! ice and snw temperatures
	128	DO jl = 1, jpl
	129	DO jk = 1, nlay_i
	130	t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
	131	END DO
	132	DO jk = 1, nlay_s
	133	t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
	134	END DO
	135	END DO
	136	!
	137	! specific temperatures for coupled runs
	138	tn_ice (:,:,:) = t_i (:,:,1,:)
	139	t1_ice (:,:,:) = t_i (:,:,1,:)
[8586]	140
[11536]	141	! heat contents
	142	e_i (:,:,:,:) = 0._wp
	143	e_s (:,:,:,:) = 0._wp
	144
	145	! general fields
	146	a_i (:,:,:) = 0._wp
	147	v_i (:,:,:) = 0._wp
	148	v_s (:,:,:) = 0._wp
	149	sv_i(:,:,:) = 0._wp
	150	oa_i(:,:,:) = 0._wp
	151	!
	152	h_i (:,:,:) = 0._wp
	153	h_s (:,:,:) = 0._wp
	154	s_i (:,:,:) = 0._wp
	155	o_i (:,:,:) = 0._wp
	156	!
	157	! melt ponds
	158	a_ip (:,:,:) = 0._wp
	159	v_ip (:,:,:) = 0._wp
	160	a_ip_frac(:,:,:) = 0._wp
	161	h_ip (:,:,:) = 0._wp
	162	!
	163	! ice velocities
	164	u_ice (:,:) = 0._wp
	165	v_ice (:,:) = 0._wp
	166	!
	167	!------------------------------------------------------------------------
	168	! 2) overwrite some of the fields with namelist parameters or netcdf file
	169	!------------------------------------------------------------------------
[8586]	170	IF( ln_iceini ) THEN
	171	! !---------------!
[11918]	172	IF( nn_iceini_file == 1 )THEN ! Read a file !
[8586]	173	! !---------------!
[11536]	174	WHERE( ff_t(:,:) >= 0._wp ) ; zswitch(:,:) = 1._wp
	175	ELSEWHERE ; zswitch(:,:) = 0._wp
	176	END WHERE
[8586]	177	!
[11536]	178	CALL fld_read( kt, 1, si ) ! input fields provided at the current time-step
[8586]	179	!
[11536]	180	! -- mandatory fields -- !
	181	zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1)
	182	zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1)
	183	zat_i_ini(:,:) = si(jp_ati)%fnow(:,:,1)
	184
	185	! -- optional fields -- !
	186	! if fields do not exist then set them to the values present in the namelist (except for snow and surface temperature)
	187	!
	188	! ice salinity
	189	IF( TRIM(si(jp_smi)%clrootname) == 'NOT USED' ) &
	190	& si(jp_smi)%fnow(:,:,1) = ( rn_smi_ini_n * zswitch + rn_smi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
	191	zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1)
	192	!
	193	! ice temperature
	194	IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) &
	195	& si(jp_tmi)%fnow(:,:,1) = ( rn_tmi_ini_n * zswitch + rn_tmi_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
	196	ztm_i_ini(:,:) = si(jp_tmi)%fnow(:,:,1)
	197	!
	198	! surface temperature => set to ice temperature if it exists
	199	IF ( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN
	200	si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
	201	ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN
	202	si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
	203	ENDIF
	204	zt_su_ini(:,:) = si(jp_tsu)%fnow(:,:,1)
	205	!
	206	! snow temperature => set to ice temperature if it exists
	207	IF ( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) == 'NOT USED' ) THEN
	208	si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
	209	ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN
	210	si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1)
	211	ENDIF
	212	ztm_s_ini(:,:) = si(jp_tms)%fnow(:,:,1)
	213	!
	214	! pond concentration
	215	IF( TRIM(si(jp_apd)%clrootname) == 'NOT USED' ) &
	216	& 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.
	217	& * si(jp_ati)%fnow(:,:,1)
	218	zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1)
	219	!
	220	! pond depth
	221	IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) &
	222	& si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1)
	223	zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1)
	224	!
	225	! change the switch for the following
	226	WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1)
	227	ELSEWHERE ; zswitch(:,:) = 0._wp
[8586]	228	END WHERE
[8813]	229	! !---------------!
[8586]	230	ELSE ! Read namelist !
	231	! !---------------!
[11536]	232	! no ice if (sst - Tfreez) >= thresold
[8813]	233	WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp
	234	ELSEWHERE ; zswitch(:,:) = tmask(:,:,1)
[8586]	235	END WHERE
[8813]	236	!
[8586]	237	! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
	238	WHERE( ff_t(:,:) >= 0._wp )
	239	zht_i_ini(:,:) = rn_hti_ini_n * zswitch(:,:)
	240	zht_s_ini(:,:) = rn_hts_ini_n * zswitch(:,:)
	241	zat_i_ini(:,:) = rn_ati_ini_n * zswitch(:,:)
	242	zsm_i_ini(:,:) = rn_smi_ini_n * zswitch(:,:)
	243	ztm_i_ini(:,:) = rn_tmi_ini_n * zswitch(:,:)
[11536]	244	zt_su_ini(:,:) = rn_tsu_ini_n * zswitch(:,:)
	245	ztm_s_ini(:,:) = rn_tms_ini_n * zswitch(:,:)
	246	zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
	247	zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:)
[8586]	248	ELSEWHERE
	249	zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:)
	250	zht_s_ini(:,:) = rn_hts_ini_s * zswitch(:,:)
	251	zat_i_ini(:,:) = rn_ati_ini_s * zswitch(:,:)
	252	zsm_i_ini(:,:) = rn_smi_ini_s * zswitch(:,:)
	253	ztm_i_ini(:,:) = rn_tmi_ini_s * zswitch(:,:)
[11536]	254	zt_su_ini(:,:) = rn_tsu_ini_s * zswitch(:,:)
	255	ztm_s_ini(:,:) = rn_tms_ini_s * zswitch(:,:)
	256	zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
	257	zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:)
[8586]	258	END WHERE
	259	!
	260	ENDIF
[11536]	261
	262	! make sure ponds = 0 if no ponds scheme
	263	IF ( .NOT.ln_pnd ) THEN
	264	zapnd_ini(:,:) = 0._wp
	265	zhpnd_ini(:,:) = 0._wp
	266	ENDIF
[8586]	267
[11536]	268	!-------------!
	269	! fill fields !
	270	!-------------!
	271	! select ice covered grid points
	272	npti = 0 ; nptidx(:) = 0
	273	DO jj = 1, jpj
	274	DO ji = 1, jpi
	275	IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
	276	npti = npti + 1
	277	nptidx(npti) = (jj - 1) * jpi + ji
	278	ENDIF
	279	END DO
	280	END DO
[8586]	281
[11536]	282	! move to 1D arrays: (jpi,jpj) -> (jpi*jpj)
	283	CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti) , zht_i_ini )
	284	CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti) , zht_s_ini )
	285	CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti) , zat_i_ini )
	286	CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini )
	287	CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini )
	288	CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti) , zt_su_ini )
	289	CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti) , zsm_i_ini )
	290	CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti) , zapnd_ini )
	291	CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti) , zhpnd_ini )
[10413]	292
[11536]	293	! allocate temporary arrays
	294	ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), &
	295	& zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) )
	296
	297	! distribute 1-cat into jpl-cat: (jpijpj) -> (jpijpj,jpl)
	298	CALL ice_var_itd( h_i_1d(1:npti) , h_s_1d(1:npti) , at_i_1d(1:npti), &
	299	& zhi_2d , zhs_2d , zai_2d , &
	300	& 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), &
	301	& zti_2d , zts_2d , ztsu_2d , zsi_2d , zaip_2d , zhip_2d )
[8586]	302
[11536]	303	! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
	304	DO jl = 1, jpl
	305	zti_3d(:,:,jl) = rt0 * tmask(:,:,1)
	306	zts_3d(:,:,jl) = rt0 * tmask(:,:,1)
	307	END DO
	308	CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d , h_i )
	309	CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d , h_s )
	310	CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d , a_i )
	311	CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d , zti_3d )
	312	CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d , zts_3d )
	313	CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d , t_su )
	314	CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d , s_i )
	315	CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d , a_ip )
	316	CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d , h_ip )
[8586]	317
[11536]	318	! deallocate temporary arrays
	319	DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
	320	& zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d )
	321
	322	! calculate extensive and intensive variables
	323	CALL ice_var_salprof ! for sz_i
	324	DO jl = 1, jpl
[8586]	325	DO jj = 1, jpj
	326	DO ji = 1, jpi
[11536]	327	v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
	328	v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
	329	sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl)
[8586]	330	END DO
	331	END DO
	332	END DO
[8813]	333	!
[11536]	334	DO jl = 1, jpl
	335	DO jk = 1, nlay_s
[8586]	336	DO jj = 1, jpj
	337	DO ji = 1, jpi
[11536]	338	t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
	339	e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
	340	& rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
[8586]	341	END DO
	342	END DO
	343	END DO
	344	END DO
[8813]	345	!
[11536]	346	DO jl = 1, jpl
	347	DO jk = 1, nlay_i
[8586]	348	DO jj = 1, jpj
	349	DO ji = 1, jpi
[11536]	350	t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl)
	351	ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
	352	e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
	353	& rhoi * ( rcpi * ( ztmelts - t_i(ji,jj,jk,jl) ) + &
	354	& rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) &
	355	& - rcp * ( ztmelts - rt0 ) )
[8586]	356	END DO
	357	END DO
	358	END DO
	359	END DO
[8933]	360
[11536]	361	! Melt ponds
	362	WHERE( a_i > epsi10 )
	363	a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
	364	ELSEWHERE
	365	a_ip_frac(:,:,:) = 0._wp
	366	END WHERE
	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	!
[11918]	471	NAMELIST/namini/ ln_iceini, nn_iceini_file, rn_thres_sst, &
[11536]	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
[11918]	498	WRITE(numout,*) ' ice initialization from a netcdf file nn_iceini_file = ', nn_iceini_file
[11536]	499	WRITE(numout,*) ' max ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst
[11927]	500	IF( ln_iceini .AND. nn_iceini_file == 0 ) THEN
[11536]	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	!
[11918]	513	IF( nn_iceini_file == 1 ) THEN ! Ice initialization using input file
[8586]	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.

Download in other formats: