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/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/ICE/iceistate.F90 @ 12644

Last change on this file since 12644 was 12644, checked in by techene, 4 years ago
stepLF: add e3 substitute and remove pe3, domqe: new routines without e3 computation
Property svn:keywords set to `Id`
File size: 26.7 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
[12583]	20	USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd
[8586]	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)
[12644]	62
[12377]	63	!! * Substitutions
	64	# include "do_loop_substitute.h90"
[8586]	65	!!----------------------------------------------------------------------
[9598]	66	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	67	!! $Id$
[10413]	68	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[8586]	69	!!----------------------------------------------------------------------
	70	CONTAINS
	71
[12377]	72	SUBROUTINE ice_istate( kt, Kbb, Kmm, Kaa )
[8586]	73	!!-------------------------------------------------------------------
	74	!! * ROUTINE ice_istate *
	75	!!
	76	!! ** Purpose : defined the sea-ice initial state
	77	!!
	78	!! ** Method : This routine will put some ice where ocean
[12583]	79	!! is at the freezing point, then fill in ice
	80	!! state variables using prescribed initial
	81	!! values in the namelist
[8586]	82	!!
	83	!! ** Steps : 1) Set initial surface and basal temperatures
	84	!! 2) Recompute or read sea ice state variables
	85	!! 3) Fill in the ice thickness distribution using gaussian
	86	!! 4) Fill in space-dependent arrays for state variables
	87	!! 5) snow-ice mass computation
	88	!! 6) store before fields
	89	!!
	90	!! ** Notes : o_i, t_su, t_s, t_i, sz_i must be filled everywhere, even
[11536]	91	!! where there is no ice
[8586]	92	!!--------------------------------------------------------------------
[12583]	93	INTEGER, INTENT(in) :: kt ! time step
[12377]	94	INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
	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(:,:) )
[12583]	119	t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1)
[11536]	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
[12583]	144
[11536]	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	! !---------------!
	172	IF( ln_iceini_file )THEN ! Read a file !
	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.
[12583]	217	& * si(jp_ati)%fnow(:,:,1)
[11536]	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
[12583]	226	WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1)
[11536]	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
[12583]	233	WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp
[8813]	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(:,:)
[12583]	246	zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.
[11536]	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
[12583]	267
[11536]	268	!-------------!
	269	! fill fields !
	270	!-------------!
	271	! select ice covered grid points
	272	npti = 0 ; nptidx(:) = 0
[12377]	273	DO_2D_11_11
	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_2D
[8586]	279
[11536]	280	! move to 1D arrays: (jpi,jpj) -> (jpi*jpj)
	281	CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti) , zht_i_ini )
	282	CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti) , zht_s_ini )
	283	CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti) , zat_i_ini )
	284	CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,1), ztm_i_ini )
	285	CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d (1:npti,1), ztm_s_ini )
	286	CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti) , zt_su_ini )
	287	CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti) , zsm_i_ini )
	288	CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti) , zapnd_ini )
	289	CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti) , zhpnd_ini )
[10413]	290
[11536]	291	! allocate temporary arrays
	292	ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), &
	293	& zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) )
[12583]	294
[11536]	295	! distribute 1-cat into jpl-cat: (jpijpj) -> (jpijpj,jpl)
	296	CALL ice_var_itd( h_i_1d(1:npti) , h_s_1d(1:npti) , at_i_1d(1:npti), &
	297	& zhi_2d , zhs_2d , zai_2d , &
	298	& 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), &
	299	& zti_2d , zts_2d , ztsu_2d , zsi_2d , zaip_2d , zhip_2d )
[8586]	300
[11536]	301	! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl)
	302	DO jl = 1, jpl
	303	zti_3d(:,:,jl) = rt0 * tmask(:,:,1)
	304	zts_3d(:,:,jl) = rt0 * tmask(:,:,1)
	305	END DO
	306	CALL tab_2d_3d( npti, nptidx(1:npti), zhi_2d , h_i )
	307	CALL tab_2d_3d( npti, nptidx(1:npti), zhs_2d , h_s )
	308	CALL tab_2d_3d( npti, nptidx(1:npti), zai_2d , a_i )
	309	CALL tab_2d_3d( npti, nptidx(1:npti), zti_2d , zti_3d )
	310	CALL tab_2d_3d( npti, nptidx(1:npti), zts_2d , zts_3d )
	311	CALL tab_2d_3d( npti, nptidx(1:npti), ztsu_2d , t_su )
	312	CALL tab_2d_3d( npti, nptidx(1:npti), zsi_2d , s_i )
	313	CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d , a_ip )
	314	CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d , h_ip )
[8586]	315
[11536]	316	! deallocate temporary arrays
	317	DEALLOCATE( zhi_2d, zhs_2d, zai_2d , &
	318	& zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d )
	319
	320	! calculate extensive and intensive variables
	321	CALL ice_var_salprof ! for sz_i
	322	DO jl = 1, jpl
[12377]	323	DO_2D_11_11
	324	v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
	325	v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
	326	sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl)
	327	END_2D
[8586]	328	END DO
[8813]	329	!
[11536]	330	DO jl = 1, jpl
[12377]	331	DO_3D_11_11( 1, nlay_s )
	332	t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
	333	e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
	334	& rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
	335	END_3D
[8586]	336	END DO
[8813]	337	!
[11536]	338	DO jl = 1, jpl
[12377]	339	DO_3D_11_11( 1, nlay_i )
[12583]	340	t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl)
[12377]	341	ztmelts = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
	342	e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
	343	& rhoi * ( rcpi * ( ztmelts - t_i(ji,jj,jk,jl) ) + &
	344	& rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) &
	345	& - rcp * ( ztmelts - rt0 ) )
	346	END_3D
[8586]	347	END DO
[8933]	348
[11536]	349	! Melt ponds
	350	WHERE( a_i > epsi10 )
	351	a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:)
	352	ELSEWHERE
	353	a_ip_frac(:,:,:) = 0._wp
	354	END WHERE
	355	v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
[12583]	356
[11536]	357	! specific temperatures for coupled runs
	358	tn_ice(:,:,:) = t_su(:,:,:)
	359	t1_ice(:,:,:) = t_i (:,:,1,:)
[8813]	360	!
[8586]	361	ENDIF ! ln_iceini
[8813]	362	!
[11536]	363	at_i(:,:) = SUM( a_i, dim=3 )
[8586]	364	!
	365	!----------------------------------------------
[11536]	366	! 3) Snow-ice mass (case ice is fully embedded)
[8586]	367	!----------------------------------------------
[9935]	368	snwice_mass (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3 ) ! snow+ice mass
[8586]	369	snwice_mass_b(:,:) = snwice_mass(:,:)
	370	!
	371	IF( ln_ice_embd ) THEN ! embedded sea-ice: deplete the initial ssh below sea-ice area
[8813]	372	!
[12377]	373	ssh(:,:,Kmm) = ssh(:,:,Kmm) - snwice_mass(:,:) * r1_rau0
	374	ssh(:,:,Kbb) = ssh(:,:,Kbb) - snwice_mass(:,:) * r1_rau0
[8813]	375	!
[12583]	376	IF( .NOT.ln_linssh ) CALL dom_vvl_zgr( Kbb, Kmm, Kaa ) ! interpolation scale factor, depth and water column
	377	! !!st
	378	! IF( .NOT.ln_linssh ) THEN
	379	! !
	380	! WHERE( ht_0(:,:) > 0 ) ; z2d(:,:) = 1._wp + ssh(:,:,Kmm)*tmask(:,:,1) / ht_0(:,:)
	381	! ELSEWHERE ; z2d(:,:) = 1._wp ; END WHERE
	382	! !
	383	! DO jk = 1,jpkm1 ! adjust initial vertical scale factors
	384	! e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * z2d(:,:)
	385	! e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm)
	386	! e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kmm)
	387	! END DO
	388	! !
	389	! ! Reconstruction of all vertical scale factors at now and before time-steps
	390	! ! =========================================================================
	391	! ! Horizontal scale factor interpolations
	392	! ! --------------------------------------
	393	! CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )
	394	! CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )
	395	! CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
	396	! CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
	397	! CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )
	398	! ! Vertical scale factor interpolations
	399	! ! ------------------------------------
	400	! CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' )
	401	! CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
	402	! CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
	403	! CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
	404	! CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
	405	! ! t- and w- points depth
	406	! ! ----------------------
	407	! !!gm not sure of that....
	408	! gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
	409	! gdepw(:,:,1,Kmm) = 0.0_wp
	410	! gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
	411	! DO jk = 2, jpk
	412	! gdept(:,:,jk,Kmm) = gdept(:,:,jk-1,Kmm) + e3w(:,:,jk ,Kmm)
	413	! gdepw(:,:,jk,Kmm) = gdepw(:,:,jk-1,Kmm) + e3t(:,:,jk-1,Kmm)
	414	! gde3w(:,:,jk) = gdept(:,:,jk ,Kmm) - ssh (:,:,Kmm)
	415	! END DO
	416	! ENDIF
[8586]	417	ENDIF
[12583]	418
[8586]	419	!------------------------------------
[11536]	420	! 4) store fields at before time-step
[8586]	421	!------------------------------------
	422	! it is only necessary for the 1st interpolation by Agrif
	423	a_i_b (:,:,:) = a_i (:,:,:)
	424	e_i_b (:,:,:,:) = e_i (:,:,:,:)
	425	v_i_b (:,:,:) = v_i (:,:,:)
	426	v_s_b (:,:,:) = v_s (:,:,:)
	427	e_s_b (:,:,:,:) = e_s (:,:,:,:)
	428	sv_i_b (:,:,:) = sv_i (:,:,:)
	429	oa_i_b (:,:,:) = oa_i (:,:,:)
	430	u_ice_b(:,:) = u_ice(:,:)
	431	v_ice_b(:,:) = v_ice(:,:)
[10527]	432	! total concentration is needed for Lupkes parameterizations
[12583]	433	at_i_b (:,:) = at_i (:,:)
[8586]	434
[8885]	435	!!clem: output of initial state should be written here but it is impossible because
	436	!! the ocean and ice are in the same file
[10425]	437	!! CALL dia_wri_state( 'output.init' )
[8813]	438	!
[8586]	439	END SUBROUTINE ice_istate
	440
[8813]	441
[8586]	442	SUBROUTINE ice_istate_init
	443	!!-------------------------------------------------------------------
	444	!! * ROUTINE ice_istate_init *
	445	!!
[12583]	446	!! ** Purpose : Definition of initial state of the ice
	447	!!
	448	!! ** Method : Read the namini namelist and check the parameter
[8813]	449	!! values called at the first timestep (nit000)
[8586]	450	!!
[8813]	451	!! ** input : Namelist namini
[8586]	452	!!
	453	!!-----------------------------------------------------------------------------
[11229]	454	INTEGER :: ios ! Local integer output status for namelist read
[8813]	455	INTEGER :: ifpr, ierror
[8586]	456	!
	457	CHARACTER(len=256) :: cn_dir ! Root directory for location of ice files
[11536]	458	TYPE(FLD_N) :: sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd
[8586]	459	TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read
	460	!
[11536]	461	NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, &
	462	& rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, &
	463	& rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, &
	464	& rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, &
	465	& rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, &
	466	& sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, cn_dir
[8586]	467	!!-----------------------------------------------------------------------------
	468	!
	469	READ ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
[11536]	470	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist' )
[8586]	471	READ ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
[11536]	472	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist' )
[8586]	473	IF(lwm) WRITE ( numoni, namini )
[8813]	474	!
[8586]	475	slf_i(jp_hti) = sn_hti ; slf_i(jp_hts) = sn_hts
[11536]	476	slf_i(jp_ati) = sn_ati ; slf_i(jp_smi) = sn_smi
	477	slf_i(jp_tmi) = sn_tmi ; slf_i(jp_tsu) = sn_tsu ; slf_i(jp_tms) = sn_tms
	478	slf_i(jp_apd) = sn_apd ; slf_i(jp_hpd) = sn_hpd
[8586]	479	!
	480	IF(lwp) THEN ! control print
	481	WRITE(numout,*)
	482	WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
	483	WRITE(numout,*) '~~~~~~~~~~~~~~~'
	484	WRITE(numout,*) ' Namelist namini:'
[11536]	485	WRITE(numout,*) ' ice initialization (T) or not (F) ln_iceini = ', ln_iceini
	486	WRITE(numout,*) ' ice initialization from a netcdf file ln_iceini_file = ', ln_iceini_file
	487	WRITE(numout,*) ' max ocean temp. above Tfreeze with initial ice rn_thres_sst = ', rn_thres_sst
	488	IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN
[12583]	489	WRITE(numout,*) ' initial snw thickness in the north-south rn_hts_ini = ', rn_hts_ini_n,rn_hts_ini_s
[11536]	490	WRITE(numout,*) ' initial ice thickness in the north-south rn_hti_ini = ', rn_hti_ini_n,rn_hti_ini_s
	491	WRITE(numout,*) ' initial ice concentr in the north-south rn_ati_ini = ', rn_ati_ini_n,rn_ati_ini_s
	492	WRITE(numout,*) ' initial ice salinity in the north-south rn_smi_ini = ', rn_smi_ini_n,rn_smi_ini_s
	493	WRITE(numout,*) ' initial surf temperat in the north-south rn_tsu_ini = ', rn_tsu_ini_n,rn_tsu_ini_s
	494	WRITE(numout,*) ' initial ice temperat in the north-south rn_tmi_ini = ', rn_tmi_ini_n,rn_tmi_ini_s
	495	WRITE(numout,*) ' initial snw temperat in the north-south rn_tms_ini = ', rn_tms_ini_n,rn_tms_ini_s
	496	WRITE(numout,*) ' initial pnd fraction in the north-south rn_apd_ini = ', rn_apd_ini_n,rn_apd_ini_s
	497	WRITE(numout,*) ' initial pnd depth in the north-south rn_hpd_ini = ', rn_hpd_ini_n,rn_hpd_ini_s
	498	ENDIF
[8586]	499	ENDIF
[8813]	500	!
[8586]	501	IF( ln_iceini_file ) THEN ! Ice initialization using input file
	502	!
	503	! set si structure
	504	ALLOCATE( si(jpfldi), STAT=ierror )
	505	IF( ierror > 0 ) THEN
[11536]	506	CALL ctl_stop( 'ice_istate_ini in iceistate: unable to allocate si structure' ) ; RETURN
[8586]	507	ENDIF
[8813]	508	!
[8586]	509	DO ifpr = 1, jpfldi
	510	ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
[11536]	511	IF( slf_i(ifpr)%ln_tint ) ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
[8586]	512	END DO
[8813]	513	!
[8586]	514	! fill si with slf_i and control print
[11536]	515	CALL fld_fill( si, slf_i, cn_dir, 'ice_istate_ini', 'initialization of sea ice fields', 'numnam_ice' )
[8813]	516	!
[8586]	517	ENDIF
[8813]	518	!
[11536]	519	IF( .NOT.ln_pnd ) THEN
	520	rn_apd_ini_n = 0. ; rn_apd_ini_s = 0.
	521	rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0.
	522	CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 when no ponds' )
	523	ENDIF
	524	!
[8586]	525	END SUBROUTINE ice_istate_init
	526
	527	#else
	528	!!----------------------------------------------------------------------
[9570]	529	!! Default option : Empty module NO SI3 sea-ice model
[8586]	530	!!----------------------------------------------------------------------
	531	#endif
	532
	533	!!======================================================================
	534	END MODULE iceistate

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

Download in other formats: