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.

icestp.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/icestp.F90 @ 8426

Last change on this file since 8426 was 8426, checked in by clem, 7 years ago
last routine names to be changed
File size: 27.9 KB

Rev	Line
[8321]	1	MODULE icestp
	2	!!======================================================================
	3	!! * MODULE icestp *
	4	!! Surface module : update the ocean surface boundary condition over ice
	5	!! & covered area using LIM sea-ice model
	6	!! Sea-Ice model : LIM-3 Sea ice model time-stepping
	7	!!=====================================================================
	8	!! History : 2.0 ! 2006-12 (M. Vancoppenolle) Original code
	9	!! 3.0 ! 2008-02 (C. Talandier) Surface module from icestp.F90
[8411]	10	!! - ! 2008-04 (G. Madec) sltyle and ice_ctl routine
[8321]	11	!! 3.3 ! 2010-11 (G. Madec) ice-ocean stress always computed at each ocean time-step
	12	!! 3.4 ! 2011-01 (A Porter) dynamical allocation
[8426]	13	!! - ! 2012-10 (C. Rousset) add ice_dia
[8321]	14	!! 3.6 ! 2014-07 (M. Vancoppenolle, G. Madec, O. Marti) revise coupled interface
	15	!! 4.0 ! 2016-06 (L. Brodeau) new unified bulk routine (based on AeroBulk)
	16	!!----------------------------------------------------------------------
	17	#if defined key_lim3
	18	!!----------------------------------------------------------------------
	19	!! 'key_lim3' : LIM 3.0 sea-ice model
	20	!!----------------------------------------------------------------------
	21	!! ice_stp : sea-ice model time-stepping and update ocean sbc over ice-covered area
	22	!!----------------------------------------------------------------------
	23	USE oce ! ocean dynamics and tracers
	24	USE dom_oce ! ocean space and time domain
	25	USE ice ! LIM-3: ice variables
[8420]	26	USE ice1D ! LIM-3: thermodynamical variables
[8321]	27	!
	28	USE sbc_oce ! Surface boundary condition: ocean fields
	29	USE sbc_ice ! Surface boundary condition: ice fields
[8404]	30	USE iceforcing ! Surface boundary condition for sea ice
[8321]	31	!
	32	USE phycst ! Define parameters for the routines
	33	USE eosbn2 ! equation of state
[8409]	34	USE icerhg ! Ice rheology
	35	USE iceadv ! Ice advection
[8422]	36	USE icethd ! Ice thermodynamics
[8409]	37	USE icerdgrft ! Ice ridging/rafting
[8414]	38	USE iceupdate ! sea surface boundary condition
[8426]	39	USE icedia ! Ice budget diagnostics
[8413]	40	USE icewri ! Ice outputs
	41	USE icerst ! Ice restarts
[8426]	42	USE icecor ! Ice corrections
[8424]	43	USE icevar ! Ice variables switch
[8411]	44	USE icectl !
[8321]	45	! MV MP 2016
	46	USE limmp
	47	! END MV MP 2016
[8426]	48	USE iceist ! LIM initial state
[8422]	49	USE icethd_sal ! LIM ice thermodynamics: salinity
[8321]	50	!
	51	USE c1d ! 1D vertical configuration
	52	USE in_out_manager ! I/O manager
	53	USE iom ! I/O manager library
	54	USE prtctl ! Print control
	55	USE lib_fortran !
	56	USE lbclnk ! lateral boundary condition - MPP link
	57	USE lib_mpp ! MPP library
	58	USE timing ! Timing
	59
	60	USE bdy_oce , ONLY: ln_bdy
[8407]	61	USE bdyice ! unstructured open boundary data
[8321]	62	# if defined key_agrif
	63	USE agrif_ice
	64	USE agrif_lim3_update
	65	USE agrif_lim3_interp
	66	# endif
	67
	68	IMPLICIT NONE
	69	PRIVATE
	70
	71	PUBLIC ice_stp ! routine called by sbcmod.F90
	72	PUBLIC ice_init ! routine called by sbcmod.F90
	73
	74	!! * Substitutions
	75	# include "vectopt_loop_substitute.h90"
	76	!!----------------------------------------------------------------------
	77	!! NEMO/OPA 4.0 , UCL NEMO Consortium (2011)
	78	!! $Id: icestp.F90 8319 2017-07-11 15:00:44Z clem $
	79	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	80	!!----------------------------------------------------------------------
	81	CONTAINS
	82
	83	SUBROUTINE ice_stp( kt, ksbc )
	84	!!---------------------------------------------------------------------
	85	!! * ROUTINE ice_stp *
	86	!!
	87	!! ** Purpose : update the ocean surface boundary condition via the
	88	!! Louvain la Neuve Sea Ice Model time stepping
	89	!!
	90	!! ** Method : ice model time stepping
	91	!! - call the ice dynamics routine
	92	!! - call the ice advection/diffusion routine
	93	!! - call the ice thermodynamics routine
	94	!! - call the routine that computes mass and
	95	!! heat fluxes at the ice/ocean interface
	96	!! - save the outputs
	97	!! - save the outputs for restart when necessary
	98	!!
	99	!! ** Action : - time evolution of the LIM sea-ice model
	100	!! - update all sbc variables below sea-ice:
	101	!! utau, vtau, taum, wndm, qns , qsr, emp , sfx
	102	!!---------------------------------------------------------------------
	103	INTEGER, INTENT(in) :: kt ! ocean time step
	104	INTEGER, INTENT(in) :: ksbc ! type of sbc flux ( 1 = user defined formulation,
	105	! 3 = bulk formulation,
	106	! 4 = Pure Coupled formulation)
	107	INTEGER :: jl ! dummy loop index
	108	!!----------------------------------------------------------------------
	109
	110	IF( nn_timing == 1 ) CALL timing_start('ice_stp')
	111
	112	!-----------------------!
	113	! --- Ice time step --- !
	114	!-----------------------!
	115	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
	116
	117	! Ice model time step
	118	kt_ice = kt
	119
	120	# if defined key_agrif
	121	IF( .NOT. Agrif_Root() ) lim_nbstep = MOD( lim_nbstep, Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) + 1
	122	# endif
	123
	124	! mean surface ocean current at ice velocity point (C-grid dynamics : U- & V-points as the ocean)
	125	u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)
	126	v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)
	127
	128	! masked sea surface freezing temperature [Kelvin] (set to rt0 over land)
	129	CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )
	130	t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
	131	!
[8404]	132	CALL ice_bef ! Store previous ice values
[8321]	133	!------------------------------------------------!
	134	! --- Dynamical coupling with the atmosphere --- !
	135	!------------------------------------------------!
[8404]	136	! it provides:
	137	! utau_ice, vtau_ice = surface ice stress [N/m2]
	138	!--------------------------------------------------
	139	CALL ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice )
	140
[8321]	141	!-------------------------------------------------------!
	142	! --- ice dynamics and transport (except in 1D case) ---!
	143	!-------------------------------------------------------!
[8404]	144	CALL ice_diag0 ! set diag of mass, heat and salt fluxes to 0
[8413]	145	CALL ice_rst_opn( kt ) ! Open Ice restart file
[8321]	146	!
	147	! --- zap this if no ice dynamics --- !
	148	IF( .NOT. lk_c1d .AND. ln_limdyn ) THEN
[8409]	149	CALL ice_rhg( kt ) ! -- rheology
	150	CALL ice_adv( kt ) ! -- advection
	151	IF( nn_limdyn == 2 .AND. nn_monocat /= 2 ) & ! -- ridging/rafting
[8426]	152	& CALL ice_rdgrft( kt )
	153	IF( nn_limdyn == 2 ) CALL ice_cor( kt , 1 ) ! -- Corrections
[8321]	154	!
	155	ENDIF
	156	! ---
[8426]	157
[8321]	158	#if defined key_agrif
	159	IF( .NOT. Agrif_Root() ) CALL agrif_interp_lim3('T')
	160	#endif
[8407]	161	IF( ln_limthd .AND. ln_bdy ) CALL bdy_ice( kt ) ! -- bdy ice thermo
[8321]	162	! previous lead fraction and ice volume for flux calculations
[8424]	163	CALL ice_var_glo2eqv ! ht_i and ht_s for ice albedo calculation
	164	CALL ice_var_agg(1) ! at_i for coupling
[8322]	165	CALL ice_bef
[8321]	166
	167	!------------------------------------------------------!
	168	! --- Thermodynamical coupling with the atmosphere --- !
	169	!------------------------------------------------------!
[8404]	170	! It provides the following fields used in sea ice model:
	171	! fr1_i0 , fr2_i0 = 1sr & 2nd fraction of qsr penetration in ice [%]
	172	! emp_oce , emp_ice = E-P over ocean and sea ice [Kg/m2/s]
	173	! sprecip = solid precipitation [Kg/m2/s]
	174	! evap_ice = sublimation [Kg/m2/s]
	175	! qsr_tot , qns_tot = solar & non solar heat flux (total) [W/m2]
	176	! qsr_ice , qns_ice = solar & non solar heat flux over ice [W/m2]
	177	! dqns_ice = non solar heat sensistivity [W/m2]
	178	! qemp_oce, qemp_ice, qprec_ice, qevap_ice = sensible heat (associated with evap & precip) [W/m2]
	179	!------------------------------------------------------------------------------------------------------
	180	CALL ice_forcing_flx( kt, ksbc )
[8321]	181
	182	!----------------------------!
	183	! --- ice thermodynamics --- !
	184	!----------------------------!
	185	! --- zap this if no ice thermo --- !
[8422]	186	IF( ln_limthd ) CALL ice_thd( kt ) ! -- Ice thermodynamics
[8321]	187
	188	! MV MP 2016
	189	IF ( ln_pnd ) CALL lim_mp( kt ) ! -- Melt ponds
	190	! END MV MP 2016
	191
[8426]	192	IF( ln_limthd ) CALL ice_cor( kt , 2 ) ! -- Corrections
[8321]	193	! ---
	194	# if defined key_agrif
	195	IF( .NOT. Agrif_Root() ) CALL agrif_update_lim3( kt )
	196	# endif
[8424]	197	CALL ice_var_glo2eqv ! necessary calls (at least for coupling)
	198	CALL ice_var_agg( 2 ) ! necessary calls (at least for coupling)
[8321]	199	!
	200	!! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work
	201	!! moreover it should only be called at the update frequency (as in agrif_lim3_update.F90)
	202	!!# if defined key_agrif
	203	!! IF( .NOT. Agrif_Root() ) CALL Agrif_ChildGrid_To_ParentGrid()
	204	!!# endif
[8414]	205	CALL ice_update_flx( kt ) ! -- Update surface ocean mass, heat and salt fluxes
[8321]	206	!!# if defined key_agrif
	207	!! IF( .NOT. Agrif_Root() ) CALL Agrif_ParentGrid_To_ChildGrid()
	208	!!# endif
[8426]	209	IF( ln_limdiahsb ) CALL ice_dia( kt ) ! -- Diagnostics and outputs
[8321]	210	!
[8413]	211	CALL ice_wri( 1 ) ! -- Ice outputs
[8321]	212	!
	213	IF( kt == nit000 .AND. ln_rstart ) &
	214	& CALL iom_close( numrir ) ! close input ice restart file
	215	!
[8413]	216	IF( lrst_ice ) CALL ice_rst_write( kt ) ! -- Ice restart file
[8321]	217	!
[8411]	218	IF( ln_limctl ) CALL ice_ctl( kt ) ! alerts in case of model crash
[8321]	219	!
	220	ENDIF ! End sea-ice time step only
	221
	222	!-------------------------!
	223	! --- Ocean time step --- !
	224	!-------------------------!
	225	! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere
	226	! using before instantaneous surf. currents
[8414]	227	IF( ln_limdyn ) CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) )
[8321]	228	!!gm remark, the ocean-ice stress is not saved in ice diag call above ..... find a solution!!!
	229	!
	230	IF( nn_timing == 1 ) CALL timing_stop('ice_stp')
	231	!
	232	END SUBROUTINE ice_stp
	233
	234
	235	SUBROUTINE ice_init
	236	!!----------------------------------------------------------------------
	237	!! * ROUTINE ice_init *
	238	!!
	239	!! ** purpose : Allocate all the dynamic arrays of the LIM-3 modules
	240	!!----------------------------------------------------------------------
	241	INTEGER :: ji, jj, ierr
	242	!!----------------------------------------------------------------------
	243	IF(lwp) WRITE(numout,*)
	244	IF(lwp) WRITE(numout,*) 'ice_init : update ocean surface boudary condition'
	245	IF(lwp) WRITE(numout,*) '~~~~~~~~ via Louvain la Neuve Ice Model (LIM-3) time stepping'
	246	!
	247	! ! Open the reference and configuration namelist files and namelist output file
	248	CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
	249	CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
	250	IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
	251	!
	252	CALL ice_run_init ! set some ice run parameters
	253	!
[8324]	254	! ! Allocate the ice arrays (sbc_ice already allocated in sbc_init)
[8321]	255	ierr = ice_alloc () ! ice variables
[8331]	256	ierr = ierr + sbc_ice_alloc () ! surface forcing
[8420]	257	ierr = ierr + ice1D_alloc () ! thermodynamics
[8409]	258	IF( ln_limdyn ) ierr = ierr + ice_rdgrft_alloc () ! ridging/rafting
[8321]	259	!
	260	IF( lk_mpp ) CALL mpp_sum( ierr )
	261	IF( ierr /= 0 ) CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays')
	262	!
[8407]	263	CALL ice_rhg_init ! set ice dynamics parameters
[8321]	264	!
	265	CALL ice_itd_init ! ice thickness distribution initialization
	266	!
[8422]	267	CALL ice_thd_init ! set ice thermodynics parameters
[8321]	268	!
[8422]	269	CALL ice_thd_sal_init ! set ice salinity parameters
[8321]	270
	271	! MV MP 2016
	272	CALL lim_mp_init ! set melt ponds parameters
	273	! END MV MP 2016
	274
[8409]	275	IF( ln_limdyn ) CALL ice_rdgrft_init ! ice thickness distribution initialization for ridging/rafting
[8321]	276	! ! Initial sea-ice state
	277	IF( .NOT. ln_rstart ) THEN ! start from rest: sea-ice deduced from sst
[8426]	278	CALL ice_ist
[8321]	279	ELSE ! start from a restart file
[8413]	280	CALL ice_rst_read
[8321]	281	ENDIF
[8424]	282	CALL ice_var_agg(2)
	283	CALL ice_var_glo2eqv
[8321]	284	!
[8414]	285	CALL ice_update_init ! ice surface boundary condition
[8321]	286	!
[8426]	287	IF( ln_limdiahsb) CALL ice_dia_init ! initialization for diags
[8321]	288	!
	289	fr_i(:,:) = at_i(:,:) ! initialisation of sea-ice fraction
	290	tn_ice(:,:,:) = t_su(:,:,:) ! initialisation of surface temp for coupled simu
	291	!
	292	DO jj = 1, jpj
	293	DO ji = 1, jpi
	294	IF( gphit(ji,jj) > 0._wp ) THEN ; rn_amax_2d(ji,jj) = rn_amax_n ! NH
	295	ELSE ; rn_amax_2d(ji,jj) = rn_amax_s ! SH
	296	ENDIF
	297	END DO
	298	END DO
	299	!
	300	END SUBROUTINE ice_init
	301
	302
	303	SUBROUTINE ice_run_init
	304	!!-------------------------------------------------------------------
	305	!! * ROUTINE ice_run_init *
	306	!!
	307	!! ** Purpose : Definition some run parameter for ice model
	308	!!
	309	!! ** Method : Read the namicerun namelist and check the parameter
	310	!! values called at the first timestep (nit000)
	311	!!
	312	!! ** input : Namelist namicerun
	313	!!-------------------------------------------------------------------
	314	INTEGER :: ios ! Local integer output status for namelist read
	315	NAMELIST/namicerun/ jpl, nlay_i, nlay_s, nn_monocat, rn_amax_n, rn_amax_s, cn_icerst_in, cn_icerst_indir, &
	316	& cn_icerst_out, cn_icerst_outdir, ln_limthd, ln_limdyn, nn_limdyn, rn_uice, rn_vice
	317	NAMELIST/namicediag/ ln_limdiachk, ln_limdiahsb, ln_limctl, iiceprt, jiceprt
	318	!!-------------------------------------------------------------------
	319	!
	320	REWIND( numnam_ice_ref ) ! Namelist namicerun in reference namelist : Parameters for ice
	321	READ ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901)
	322	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp )
	323
	324	REWIND( numnam_ice_cfg ) ! Namelist namicerun in configuration namelist : Parameters for ice
	325	READ ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 )
	326	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp )
	327	IF(lwm) WRITE ( numoni, namicerun )
	328	!
	329	REWIND( numnam_ice_ref ) ! Namelist namicediag in reference namelist : Parameters for ice
	330	READ ( numnam_ice_ref, namicediag, IOSTAT = ios, ERR = 903)
	331	903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in reference namelist', lwp )
	332
	333	REWIND( numnam_ice_cfg ) ! Namelist namicediag in configuration namelist : Parameters for ice
	334	READ ( numnam_ice_cfg, namicediag, IOSTAT = ios, ERR = 904 )
	335	904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in configuration namelist', lwp )
	336	IF(lwm) WRITE ( numoni, namicediag )
	337	!
	338	IF(lwp) THEN ! control print
	339	WRITE(numout,*)
	340	WRITE(numout,*) 'ice_run_init : ice share parameters for dynamics/advection/thermo of sea-ice'
	341	WRITE(numout,*) ' ~~~~~~'
	342	WRITE(numout,*) ' number of ice categories jpl = ', jpl
	343	WRITE(numout,*) ' number of ice layers nlay_i = ', nlay_i
	344	WRITE(numout,*) ' number of snow layers nlay_s = ', nlay_s
	345	WRITE(numout,*) ' virtual ITD mono-category param (1-4) or not (0) nn_monocat = ', nn_monocat
	346	WRITE(numout,*) ' maximum ice concentration for NH = ', rn_amax_n
	347	WRITE(numout,*) ' maximum ice concentration for SH = ', rn_amax_s
	348	WRITE(numout,*) ' Ice thermodynamics (T) or not (F) ln_limthd = ', ln_limthd
	349	WRITE(numout,*) ' Ice dynamics (T) or not (F) ln_limdyn = ', ln_limdyn
	350	WRITE(numout,*) ' (ln_limdyn=T) Ice dynamics switch nn_limdyn = ', nn_limdyn
	351	WRITE(numout,*) ' 2: total'
	352	WRITE(numout,*) ' 1: advection only (no diffusion, no ridging/rafting)'
	353	WRITE(numout,*) ' 0: advection only (as 1 + prescribed velocity, bypass rheology)'
	354	WRITE(numout,*) ' (ln_limdyn=T) prescribed u-vel (case nn_limdyn=0) = ', rn_uice
	355	WRITE(numout,*) ' (ln_limdyn=T) prescribed v-vel (case nn_limdyn=0) = ', rn_vice
	356	WRITE(numout,*)
	357	WRITE(numout,*) '...and ice diagnostics'
	358	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~'
	359	WRITE(numout,*) ' Diagnose online heat/mass/salt budget ln_limdiachk = ', ln_limdiachk
	360	WRITE(numout,*) ' Output heat/mass/salt budget ln_limdiahsb = ', ln_limdiahsb
	361	WRITE(numout,*) ' control prints in ocean.out for (i,j)=(iiceprt,jiceprt) = ', ln_limctl
	362	WRITE(numout,*) ' i-index for control prints (ln_limctl=true) = ', iiceprt
	363	WRITE(numout,*) ' j-index for control prints (ln_limctl=true) = ', jiceprt
	364	ENDIF
	365	!
	366	IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN
	367	nn_monocat = 0
	368	IF(lwp) WRITE(numout,*)
	369	IF(lwp) WRITE(numout,*) ' nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
	370	ENDIF
	371	IF ( ( jpl == 1 ) .AND. ( nn_monocat == 0 ) ) THEN
	372	CALL ctl_stop( 'STOP', 'ice_run_init : if jpl=1 then nn_monocat should be between 1 and 4' )
	373	ENDIF
	374	!
	375	! sea-ice timestep and inverse
	376	rdt_ice = REAL(nn_fsbc) * rdt
	377	r1_rdtice = 1._wp / rdt_ice
	378
	379	! inverse of nlay_i and nlay_s
	380	r1_nlay_i = 1._wp / REAL( nlay_i, wp )
	381	r1_nlay_s = 1._wp / REAL( nlay_s, wp )
	382	!
	383	IF( lwp .AND. ln_bdy .AND. ln_limdiachk ) &
	384	& CALL ctl_warn('online conservation check activated but it does not work with BDY')
	385	!
	386	IF( lwp ) WRITE(numout,*) ' ice timestep rdt_ice = ', rdt_ice
	387	!
	388	END SUBROUTINE ice_run_init
	389
	390
	391	SUBROUTINE ice_itd_init
	392	!!------------------------------------------------------------------
	393	!! * ROUTINE ice_itd_init *
	394	!!
	395	!! ** Purpose : Initializes the ice thickness distribution
	396	!! ** Method : ...
	397	!! ** input : Namelist namiceitd
	398	!!-------------------------------------------------------------------
	399	INTEGER :: ios ! Local integer output status for namelist read
	400	NAMELIST/namiceitd/ rn_himean
	401	!
	402	INTEGER :: jl ! dummy loop index
	403	REAL(wp) :: zc1, zc2, zc3, zx1 ! local scalars
	404	REAL(wp) :: zhmax, znum, zden, zalpha !
	405	!!------------------------------------------------------------------
	406	!
	407	REWIND( numnam_ice_ref ) ! Namelist namiceitd in reference namelist : Parameters for ice
	408	READ ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 905)
	409	905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp )
	410
	411	REWIND( numnam_ice_cfg ) ! Namelist namiceitd in configuration namelist : Parameters for ice
	412	READ ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 906 )
	413	906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp )
	414	IF(lwm) WRITE ( numoni, namiceitd )
	415	!
	416	IF(lwp) THEN ! control print
	417	WRITE(numout,*)
	418	WRITE(numout,*) 'ice_itd_init : Initialization of ice cat distribution '
	419	WRITE(numout,*) '~~~~~~~~~~~~'
	420	WRITE(numout,*) ' mean ice thickness in the domain rn_himean = ', rn_himean
	421	ENDIF
	422	!
	423	!----------------------------------
	424	!- Thickness categories boundaries
	425	!----------------------------------
	426	!
	427	!== h^(-alpha) function ==!
	428	zalpha = 0.05_wp
	429	zhmax = 3._wp * rn_himean
	430	DO jl = 1, jpl
	431	znum = jpl * ( zhmax+1 )**zalpha
	432	zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
	433	hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
	434	END DO
	435	!
	436	DO jl = 1, jpl ! mean thickness by category
	437	hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
	438	END DO
	439	!
	440	hi_max(jpl) = 99._wp ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
	441	!
	442	IF(lwp) WRITE(numout,*)
	443	IF(lwp) WRITE(numout,*) ' Thickness category boundaries '
	444	IF(lwp) WRITE(numout,*) ' hi_max ', hi_max(0:jpl)
	445	!
	446	END SUBROUTINE ice_itd_init
	447
	448	SUBROUTINE ice_bef
	449	!!----------------------------------------------------------------------
	450	!! * ROUTINE ice_bef *
	451	!!
	452	!! ** purpose : store ice variables at "before" time step
	453	!!----------------------------------------------------------------------
[8360]	454	INTEGER :: ji, jj, jl ! dummy loop index
	455	!!----------------------------------------------------------------------
[8321]	456	!
[8360]	457	DO jl = 1, jpl
[8404]	458	DO jj = 2, jpjm1
	459	DO ji = 2, jpim1
[8360]	460	a_i_b (ji,jj,jl) = a_i (ji,jj,jl) ! ice area
	461	v_i_b (ji,jj,jl) = v_i (ji,jj,jl) ! ice volume
	462	v_s_b (ji,jj,jl) = v_s (ji,jj,jl) ! snow volume
	463	smv_i_b(ji,jj,jl) = smv_i(ji,jj,jl) ! salt content
	464	oa_i_b (ji,jj,jl) = oa_i (ji,jj,jl) ! areal age content
	465	e_s_b (ji,jj,:,jl) = e_s (ji,jj,:,jl) ! snow thermal energy
	466	e_i_b (ji,jj,:,jl) = e_i (ji,jj,:,jl) ! ice thermal energy
	467	! ! ice thickness
	468	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi20 ) ) ! 0 if no ice and 1 if yes
	469	ht_i_b(ji,jj,jl) = v_i_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch
	470	ht_s_b(ji,jj,jl) = v_s_b (ji,jj,jl) / MAX( a_i_b(ji,jj,jl) , epsi20 ) * rswitch
	471	END DO
[8404]	472	END DO
[8360]	473	END DO
[8404]	474	CALL lbc_lnk_multi( a_i_b, 'T', 1., v_i_b , 'T', 1., v_s_b , 'T', 1., smv_i_b, 'T', 1., &
	475	& oa_i_b, 'T', 1., ht_i_b, 'T', 1., ht_s_b, 'T', 1. )
	476	CALL lbc_lnk( e_i_b, 'T', 1. )
	477	CALL lbc_lnk( e_s_b, 'T', 1. )
[8321]	478
[8360]	479	! ice velocities & total concentration
[8404]	480	DO jj = 2, jpjm1
	481	DO ji = 2, jpim1
[8360]	482	at_i_b(ji,jj) = SUM( a_i_b(ji,jj,:) )
	483	u_ice_b(ji,jj) = u_ice(ji,jj)
	484	v_ice_b(ji,jj) = v_ice(ji,jj)
	485	END DO
	486	END DO
[8404]	487	CALL lbc_lnk_multi( at_i_b, 'T', 1., u_ice_b , 'U', -1., v_ice_b , 'V', -1. )
[8360]	488
[8321]	489	END SUBROUTINE ice_bef
	490
	491
	492	SUBROUTINE ice_diag0
	493	!!----------------------------------------------------------------------
	494	!! * ROUTINE ice_diag0 *
	495	!!
	496	!! ** purpose : set ice-ocean and ice-atm. fluxes to zeros at the beggining
	497	!! of the time step
	498	!!----------------------------------------------------------------------
[8360]	499	INTEGER :: ji, jj ! dummy loop index
	500	!!----------------------------------------------------------------------
	501	DO jj = 1, jpj
	502	DO ji = 1, jpi
	503	sfx (ji,jj) = 0._wp ;
	504	sfx_bri(ji,jj) = 0._wp ; sfx_lam(ji,jj) = 0._wp
	505	sfx_sni(ji,jj) = 0._wp ; sfx_opw(ji,jj) = 0._wp
	506	sfx_bog(ji,jj) = 0._wp ; sfx_dyn(ji,jj) = 0._wp
	507	sfx_bom(ji,jj) = 0._wp ; sfx_sum(ji,jj) = 0._wp
	508	sfx_res(ji,jj) = 0._wp ; sfx_sub(ji,jj) = 0._wp
	509	!
	510	wfx_snw(ji,jj) = 0._wp ; wfx_ice(ji,jj) = 0._wp
	511	wfx_sni(ji,jj) = 0._wp ; wfx_opw(ji,jj) = 0._wp
	512	wfx_bog(ji,jj) = 0._wp ; wfx_dyn(ji,jj) = 0._wp
	513	wfx_bom(ji,jj) = 0._wp ; wfx_sum(ji,jj) = 0._wp
	514	wfx_res(ji,jj) = 0._wp ; wfx_sub(ji,jj) = 0._wp
	515	wfx_spr(ji,jj) = 0._wp ; wfx_lam(ji,jj) = 0._wp
	516	wfx_snw_dyn(ji,jj) = 0._wp ; wfx_snw_sum(ji,jj) = 0._wp
	517	wfx_snw_sub(ji,jj) = 0._wp ; wfx_ice_sub(ji,jj) = 0._wp
	518	wfx_snw_sni(ji,jj) = 0._wp
	519	! MV MP 2016
	520	wfx_pnd(ji,jj) = 0._wp
	521	! END MV MP 2016
	522
	523	hfx_thd(ji,jj) = 0._wp ;
	524	hfx_snw(ji,jj) = 0._wp ; hfx_opw(ji,jj) = 0._wp
	525	hfx_bog(ji,jj) = 0._wp ; hfx_dyn(ji,jj) = 0._wp
	526	hfx_bom(ji,jj) = 0._wp ; hfx_sum(ji,jj) = 0._wp
	527	hfx_res(ji,jj) = 0._wp ; hfx_sub(ji,jj) = 0._wp
	528	hfx_spr(ji,jj) = 0._wp ; hfx_dif(ji,jj) = 0._wp
	529	hfx_err(ji,jj) = 0._wp ; hfx_err_rem(ji,jj) = 0._wp
	530	hfx_err_dif(ji,jj) = 0._wp
	531	wfx_err_sub(ji,jj) = 0._wp
	532	!
	533	afx_tot(ji,jj) = 0._wp ;
	534	afx_dyn(ji,jj) = 0._wp ; afx_thd(ji,jj) = 0._wp
	535	!
	536	diag_heat(ji,jj) = 0._wp ; diag_smvi(ji,jj) = 0._wp
	537	diag_vice(ji,jj) = 0._wp ; diag_vsnw(ji,jj) = 0._wp
	538
	539	! SIMIP diagnostics
	540	diag_dms_dyn(ji,jj) = 0._wp ; diag_dmi_dyn(ji,jj) = 0._wp
	541	diag_fc_bo(ji,jj) = 0._wp ; diag_fc_su(ji,jj) = 0._wp
	542
	543	tau_icebfr(ji,jj) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
	544	END DO
	545	END DO
[8321]	546
	547	END SUBROUTINE ice_diag0
	548
	549	#else
	550	!!----------------------------------------------------------------------
	551	!! Default option Dummy module NO LIM 3.0 sea-ice model
	552	!!----------------------------------------------------------------------
	553	CONTAINS
	554	!
	555	SUBROUTINE ice_stp ( kt, ksbc ) ! Dummy routine
	556	INTEGER, INTENT(in) :: kt, ksbc
	557	WRITE(,) 'ice_stp: You should not have seen this print! error?', kt
	558	END SUBROUTINE ice_stp
	559	!
	560	SUBROUTINE ice_init ! Dummy routine
	561	END SUBROUTINE ice_init
	562	#endif
	563
	564	!!======================================================================
	565	END MODULE icestp

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

Download in other formats: