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.

sbcice_lim.F90 in branches/UKMO/dev_r8864_restart_date/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/UKMO/dev_r8864_restart_date/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90 @ 9235

Last change on this file since 9235 was 9235, checked in by davestorkey, 6 years ago
UKMO/dev_r8864_restart_date : clear SVN keywords.
File size: 34.1 KB

Rev	Line
[888]	1	MODULE sbcice_lim
	2	!!======================================================================
	3	!! * MODULE sbcice_lim *
[921]	4	!! Surface module : update the ocean surface boundary condition over ice
	5	!! & covered area using LIM sea-ice model
[2715]	6	!! Sea-Ice model : LIM-3 Sea ice model time-stepping
[2528]	7	!!=====================================================================
	8	!! History : 2.0 ! 2006-12 (M. Vancoppenolle) Original code
	9	!! 3.0 ! 2008-02 (C. Talandier) Surface module from icestp.F90
	10	!! - ! 2008-04 (G. Madec) sltyle and lim_ctl routine
	11	!! 3.3 ! 2010-11 (G. Madec) ice-ocean stress always computed at each ocean time-step
[3625]	12	!! 3.4 ! 2011-01 (A Porter) dynamical allocation
[4161]	13	!! - ! 2012-10 (C. Rousset) add lim_diahsb
[4990]	14	!! 3.6 ! 2014-07 (M. Vancoppenolle, G. Madec, O. Marti) revise coupled interface
[7646]	15	!! 4.0 ! 2016-06 (L. Brodeau) new unified bulk routine (based on AeroBulk)
[888]	16	!!----------------------------------------------------------------------
	17	#if defined key_lim3
	18	!!----------------------------------------------------------------------
	19	!! 'key_lim3' : LIM 3.0 sea-ice model
	20	!!----------------------------------------------------------------------
[921]	21	!! sbc_ice_lim : sea-ice model time-stepping and update ocean sbc over ice-covered area
[888]	22	!!----------------------------------------------------------------------
	23	USE oce ! ocean dynamics and tracers
	24	USE dom_oce ! ocean space and time domain
[2528]	25	USE ice ! LIM-3: ice variables
[5123]	26	USE thd_ice ! LIM-3: thermodynamical variables
[6140]	27	!
[888]	28	USE sbc_oce ! Surface boundary condition: ocean fields
	29	USE sbc_ice ! Surface boundary condition: ice fields
[7646]	30	USE usrdef_sbc ! user defined: surface boundary condition
	31	USE sbcblk ! Surface boundary condition: bulk
[4161]	32	USE sbccpl ! Surface boundary condition: coupled interface
[2528]	33	USE albedo ! ocean & ice albedo
[6140]	34	!
[888]	35	USE phycst ! Define parameters for the routines
	36	USE eosbn2 ! equation of state
	37	USE limdyn ! Ice dynamics
	38	USE limtrp ! Ice transport
[5429]	39	USE limhdf ! Ice horizontal diffusion
[888]	40	USE limthd ! Ice thermodynamics
	41	USE limitd_me ! Mechanics on ice thickness distribution
	42	USE limsbc ! sea surface boundary condition
[4161]	43	USE limdiahsb ! Ice budget diagnostics
[888]	44	USE limwri ! Ice outputs
	45	USE limrst ! Ice restarts
[5123]	46	USE limupdate1 ! update of global variables
	47	USE limupdate2 ! update of global variables
[888]	48	USE limvar ! Ice variables switch
[7646]	49	USE limctl !
[5123]	50	USE limistate ! LIM initial state
	51	USE limthd_sal ! LIM ice thermodynamics: salinity
[6140]	52	!
[2528]	53	USE c1d ! 1D vertical configuration
[6140]	54	USE in_out_manager ! I/O manager
	55	USE iom ! I/O manager library
	56	USE prtctl ! Print control
[7646]	57	USE lib_fortran !
[2528]	58	USE lbclnk ! lateral boundary condition - MPP link
[2715]	59	USE lib_mpp ! MPP library
[3294]	60	USE wrk_nemo ! work arrays
[4161]	61	USE timing ! Timing
[888]	62
[7646]	63	USE bdy_oce , ONLY: ln_bdy
	64	USE bdyice_lim ! unstructured open boundary data (bdy_ice_lim routine)
	65	# if defined key_agrif
	66	USE agrif_ice
	67	USE agrif_lim3_update
	68	USE agrif_lim3_interp
	69	# endif
[4161]	70
[888]	71	IMPLICIT NONE
	72	PRIVATE
	73
	74	PUBLIC sbc_ice_lim ! routine called by sbcmod.F90
[5123]	75	PUBLIC sbc_lim_init ! routine called by sbcmod.F90
[7646]	76
[888]	77	!! * Substitutions
	78	# include "vectopt_loop_substitute.h90"
	79	!!----------------------------------------------------------------------
[2715]	80	!! NEMO/OPA 4.0 , UCL NEMO Consortium (2011)
[1146]	81	!! $Id$
[2528]	82	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[888]	83	!!----------------------------------------------------------------------
	84	CONTAINS
	85
[7646]	86	SUBROUTINE sbc_ice_lim( kt, ksbc )
[888]	87	!!---------------------------------------------------------------------
	88	!! * ROUTINE sbc_ice_lim *
	89	!!
[7646]	90	!! ** Purpose : update the ocean surface boundary condition via the
	91	!! Louvain la Neuve Sea Ice Model time stepping
	92	!!
[888]	93	!! ** Method : ice model time stepping
[7646]	94	!! - call the ice dynamics routine
	95	!! - call the ice advection/diffusion routine
	96	!! - call the ice thermodynamics routine
	97	!! - call the routine that computes mass and
[888]	98	!! heat fluxes at the ice/ocean interface
[7646]	99	!! - save the outputs
[888]	100	!! - save the outputs for restart when necessary
	101	!!
	102	!! ** Action : - time evolution of the LIM sea-ice model
	103	!! - update all sbc variables below sea-ice:
[7646]	104	!! utau, vtau, taum, wndm, qns , qsr, emp , sfx
[888]	105	!!---------------------------------------------------------------------
	106	INTEGER, INTENT(in) :: kt ! ocean time step
[7646]	107	INTEGER, INTENT(in) :: ksbc ! type of sbc flux ( 1 = user defined formulation,
	108	! 3 = bulk formulation,
	109	! 4 = Pure Coupled formulation)
[888]	110	!!
[7753]	111	INTEGER :: jl ! dummy loop index
[4990]	112	REAL(wp), POINTER, DIMENSION(:,:,:) :: zalb_os, zalb_cs ! ice albedo under overcast/clear sky
[5407]	113	REAL(wp), POINTER, DIMENSION(:,: ) :: zutau_ice, zvtau_ice
[888]	114	!!----------------------------------------------------------------------
	115
[7646]	116	IF( nn_timing == 1 ) CALL timing_start('sbc_ice_lim')
[4161]	117
[5410]	118	!-----------------------!
	119	! --- Ice time step --- !
	120	!-----------------------!
	121	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
[5407]	122
[7646]	123	# if defined key_agrif
	124	IF( .NOT. Agrif_Root() ) lim_nbstep = MOD( lim_nbstep, Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) + 1
	125	# endif
	126
[5410]	127	! mean surface ocean current at ice velocity point (C-grid dynamics : U- & V-points as the ocean)
[7753]	128	u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)
	129	v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)
[7646]	130
[5123]	131	! masked sea surface freezing temperature [Kelvin] (set to rt0 over land)
[5541]	132	CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )
[7753]	133	t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
[7646]	134
[5123]	135	! Mask sea ice surface temperature (set to rt0 over land)
[1109]	136	DO jl = 1, jpl
[7753]	137	t_su(:,:,jl) = t_su(:,:,jl) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) )
[7646]	138	END DO
[888]	139	!
[5410]	140	!------------------------------------------------!
[7646]	141	! --- Dynamical coupling with the atmosphere --- !
	142	!------------------------------------------------!
[5410]	143	! It provides the following fields:
	144	! utau_ice, vtau_ice : surface ice stress (U- & V-points) [N/m2]
	145	!-----------------------------------------------------------------
[7646]	146	CALL sbc_lim_bef ! Store previous ice values
	147	SELECT CASE( ksbc )
	148	CASE( jp_usr ) ; CALL usrdef_sbc_ice_tau( kt ) ! user defined formulation
	149	CASE( jp_blk ) ; CALL blk_ice_tau ! Bulk formulation
	150	CASE( jp_purecpl ) ; CALL sbc_cpl_ice_tau( utau_ice , vtau_ice ) ! Coupled formulation
[888]	151	END SELECT
[7646]	152
	153	IF( ln_mixcpl) THEN ! Case of a mixed Bulk/Coupled formulation
[5407]	154	CALL wrk_alloc( jpi,jpj , zutau_ice, zvtau_ice)
[7646]	155	CALL sbc_cpl_ice_tau( zutau_ice , zvtau_ice )
[7753]	156	utau_ice(:,:) = utau_ice(:,:) * xcplmask(:,:,0) + zutau_ice(:,:) * ( 1. - xcplmask(:,:,0) )
	157	vtau_ice(:,:) = vtau_ice(:,:) * xcplmask(:,:,0) + zvtau_ice(:,:) * ( 1. - xcplmask(:,:,0) )
[5407]	158	CALL wrk_dealloc( jpi,jpj , zutau_ice, zvtau_ice)
	159	ENDIF
	160
[5410]	161	!-------------------------------------------------------!
	162	! --- ice dynamics and transport (except in 1D case) ---!
	163	!-------------------------------------------------------!
	164	numit = numit + nn_fsbc ! Ice model time step
[921]	165	!
[7646]	166	CALL sbc_lim_diag0 ! set diag of mass, heat and salt fluxes to 0
	167	CALL lim_rst_opn( kt ) ! Open Ice restart file
	168	!
	169	! --- zap this if no ice dynamics --- !
	170	IF( .NOT. lk_c1d .AND. ln_limdyn ) THEN
[5410]	171	!
[7646]	172	IF( nn_limdyn /= 0 ) THEN ! -- Ice dynamics
	173	CALL lim_dyn( kt ) ! rheology
	174	ELSE
[7753]	175	u_ice(:,:) = rn_uice * umask(:,:,1) ! or prescribed velocity
	176	v_ice(:,:) = rn_vice * vmask(:,:,1)
[7646]	177	ENDIF
	178	CALL lim_trp( kt ) ! -- Ice transport (Advection/diffusion)
	179	IF( nn_limdyn == 2 .AND. nn_monocat /= 2 ) & ! -- Mechanical redistribution (ridging/rafting)
	180	& CALL lim_itd_me
	181	IF( nn_limdyn == 2 ) CALL lim_update1( kt ) ! -- Corrections
[5410]	182	!
[7646]	183	ENDIF
	184
	185	! ---
	186	#if defined key_agrif
	187	IF( .NOT. Agrif_Root() ) CALL agrif_interp_lim3('T')
[4333]	188	#endif
[7761]	189	IF( ln_limthd .AND. ln_bdy ) CALL bdy_ice_lim( kt ) ! -- bdy ice thermo
[5123]	190	! previous lead fraction and ice volume for flux calculations
[7646]	191	CALL sbc_lim_bef
	192	CALL lim_var_glo2eqv ! ht_i and ht_s for ice albedo calculation
	193	CALL lim_var_agg(1) ! at_i for coupling (via pfrld)
	194	!
[7753]	195	pfrld(:,:) = 1._wp - at_i(:,:)
	196	phicif(:,:) = vt_i(:,:)
[7646]	197
[5410]	198	!------------------------------------------------------!
[7646]	199	! --- Thermodynamical coupling with the atmosphere --- !
	200	!------------------------------------------------------!
[5410]	201	! It provides the following fields:
	202	! qsr_ice , qns_ice : solar & non solar heat flux over ice (T-point) [W/m2]
	203	! qla_ice : latent heat flux over ice (T-point) [W/m2]
	204	! dqns_ice, dqla_ice : non solar & latent heat sensistivity (T-point) [W/m2]
	205	! tprecip , sprecip : total & solid precipitation (T-point) [Kg/m2/s]
	206	! fr1_i0 , fr2_i0 : 1sr & 2nd fraction of qsr penetration in ice [%]
	207	!----------------------------------------------------------------------------------------
[6416]	208	CALL wrk_alloc( jpi,jpj,jpl, zalb_os, zalb_cs )
[7646]	209
	210	CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
	211	SELECT CASE( ksbc )
	212	CASE( jp_usr ) ; CALL usrdef_sbc_ice_flx( kt ) ! user defined formulation
	213	CASE( jp_blk ) ! bulk formulation
	214	! albedo depends on cloud fraction because of non-linear spectral effects
[7753]	215	alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
[7646]	216	CALL blk_ice_flx( t_su, alb_ice )
	217	IF( ln_mixcpl ) CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
	218	IF( nn_limflx /= 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
	219	CASE ( jp_purecpl )
	220	! albedo depends on cloud fraction because of non-linear spectral effects
[7753]	221	alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
[7646]	222	CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
	223	IF( nn_limflx == 2 ) CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
	224	END SELECT
[5410]	225
[6416]	226	CALL wrk_dealloc( jpi,jpj,jpl, zalb_os, zalb_cs )
[5407]	227
[5410]	228	!----------------------------!
	229	! --- ice thermodynamics --- !
	230	!----------------------------!
[7646]	231	! --- zap this if no ice thermo --- !
	232	IF( ln_limthd ) CALL lim_thd( kt ) ! -- Ice thermodynamics
	233	IF( ln_limthd ) CALL lim_update2( kt ) ! -- Corrections
	234	! ---
	235	# if defined key_agrif
	236	IF( .NOT. Agrif_Root() ) CALL agrif_update_lim3( kt )
	237	# endif
	238	CALL lim_var_glo2eqv ! necessary calls (at least for coupling)
	239	CALL lim_var_agg( 2 ) ! necessary calls (at least for coupling)
	240	!
[7761]	241	!! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work
	242	!! moreover it should only be called at the update frequency (as in agrif_lim3_update.F90)
	243	!!# if defined key_agrif
	244	!! IF( .NOT. Agrif_Root() ) CALL Agrif_ChildGrid_To_ParentGrid()
	245	!!# endif
[7646]	246	CALL lim_sbc_flx( kt ) ! -- Update surface ocean mass, heat and salt fluxes
[7761]	247	!!# if defined key_agrif
	248	!! IF( .NOT. Agrif_Root() ) CALL Agrif_ParentGrid_To_ChildGrid()
	249	!!# endif
[7646]	250	IF( ln_limdiahsb ) CALL lim_diahsb( kt ) ! -- Diagnostics and outputs
[921]	251	!
[7646]	252	CALL lim_wri( 1 ) ! -- Ice outputs
[888]	253	!
[4205]	254	IF( kt == nit000 .AND. ln_rstart ) &
[7646]	255	& CALL iom_close( numrir ) ! close input ice restart file
[4205]	256	!
[7646]	257	IF( lrst_ice ) CALL lim_rst_write( kt ) ! -- Ice restart file
[921]	258	!
[7646]	259	IF( ln_limctl ) CALL lim_ctl( kt ) ! alerts in case of model crash
[921]	260	!
[5123]	261	ENDIF ! End sea-ice time step only
[888]	262
[5410]	263	!-------------------------!
	264	! --- Ocean time step --- !
	265	!-------------------------!
	266	! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere
[7646]	267	! using before instantaneous surf. currents
	268	IF( ln_limdyn ) CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )
[2528]	269	!!gm remark, the ocean-ice stress is not saved in ice diag call above ..... find a solution!!!
[2715]	270	!
[5410]	271	IF( nn_timing == 1 ) CALL timing_stop('sbc_ice_lim')
[4161]	272	!
[921]	273	END SUBROUTINE sbc_ice_lim
[5123]	274
[7646]	275
[5123]	276	SUBROUTINE sbc_lim_init
	277	!!----------------------------------------------------------------------
	278	!! * ROUTINE sbc_lim_init *
	279	!!
	280	!! ** purpose : Allocate all the dynamic arrays of the LIM-3 modules
	281	!!----------------------------------------------------------------------
[7753]	282	INTEGER :: ji, jj, ierr
[5123]	283	!!----------------------------------------------------------------------
	284	IF(lwp) WRITE(numout,*)
[7646]	285	IF(lwp) WRITE(numout,*) 'sbc_lim_init : update ocean surface boudary condition'
[5123]	286	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ via Louvain la Neuve Ice Model (LIM-3) time stepping'
	287	!
[7646]	288	! ! Open the reference and configuration namelist files and namelist output file
	289	CALL ctl_opn( numnam_ice_ref, 'namelist_ice_ref', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
[5123]	290	CALL ctl_opn( numnam_ice_cfg, 'namelist_ice_cfg', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
	291	IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
[6140]	292	!
[7646]	293	CALL lim_run_init ! set some ice run parameters
[5123]	294	!
	295	! ! Allocate the ice arrays
	296	ierr = ice_alloc () ! ice variables
	297	ierr = ierr + sbc_ice_alloc () ! surface forcing
	298	ierr = ierr + thd_ice_alloc () ! thermodynamics
[7646]	299	IF( ln_limdyn ) ierr = ierr + lim_itd_me_alloc () ! ice thickness distribution - mechanics
[5123]	300	!
	301	IF( lk_mpp ) CALL mpp_sum( ierr )
	302	IF( ierr /= 0 ) CALL ctl_stop('STOP', 'sbc_lim_init : unable to allocate ice arrays')
	303	!
[7646]	304	CALL lim_dyn_init ! set ice dynamics parameters
[5123]	305	!
	306	CALL lim_itd_init ! ice thickness distribution initialization
	307	!
[5429]	308	CALL lim_hdf_init ! set ice horizontal diffusion computation parameters
	309	!
[5123]	310	CALL lim_thd_init ! set ice thermodynics parameters
	311	!
	312	CALL lim_thd_sal_init ! set ice salinity parameters
	313	!
[7646]	314	IF( ln_limdyn ) CALL lim_itd_me_init ! ice thickness distribution initialization for mecanical deformation
[5123]	315	! ! Initial sea-ice state
	316	IF( .NOT. ln_rstart ) THEN ! start from rest: sea-ice deduced from sst
	317	numit = 0
	318	numit = nit000 - 1
	319	CALL lim_istate
	320	ELSE ! start from a restart file
	321	CALL lim_rst_read
	322	numit = nit000 - 1
	323	ENDIF
[7646]	324	CALL lim_var_agg(2)
[5123]	325	CALL lim_var_glo2eqv
	326	!
[7646]	327	CALL lim_sbc_init ! ice surface boundary condition
[5123]	328	!
[7646]	329	IF( ln_limdiahsb) CALL lim_diahsb_init ! initialization for diags
	330	!
[7753]	331	fr_i(:,:) = at_i(:,:) ! initialisation of sea-ice fraction
	332	tn_ice(:,:,:) = t_su(:,:,:) ! initialisation of surface temp for coupled simu
[5123]	333	!
[6403]	334	DO jj = 1, jpj
	335	DO ji = 1, jpi
	336	IF( gphit(ji,jj) > 0._wp ) THEN ; rn_amax_2d(ji,jj) = rn_amax_n ! NH
	337	ELSE ; rn_amax_2d(ji,jj) = rn_amax_s ! SH
	338	ENDIF
[7646]	339	END DO
	340	END DO
[6403]	341	!
[7646]	342	nstart = numit + nn_fsbc
	343	nitrun = nitend - nit000 + 1
	344	nlast = numit + nitrun
[5123]	345	!
	346	IF( nstock == 0 ) nstock = nlast + 1
	347	!
	348	END SUBROUTINE sbc_lim_init
	349
	350
[7646]	351	SUBROUTINE lim_run_init
[5123]	352	!!-------------------------------------------------------------------
[7646]	353	!! * ROUTINE lim_run_init *
	354	!!
[5123]	355	!! ** Purpose : Definition some run parameter for ice model
	356	!!
[7646]	357	!! ** Method : Read the namicerun namelist and check the parameter
[5123]	358	!! values called at the first timestep (nit000)
	359	!!
	360	!! ** input : Namelist namicerun
	361	!!-------------------------------------------------------------------
	362	INTEGER :: ios ! Local integer output status for namelist read
[7646]	363	NAMELIST/namicerun/ jpl, nlay_i, nlay_s, rn_amax_n, rn_amax_s, cn_icerst_in, cn_icerst_indir, &
	364	& cn_icerst_out, cn_icerst_outdir, ln_limthd, ln_limdyn, nn_limdyn, rn_uice, rn_vice
	365	NAMELIST/namicediag/ ln_limdiachk, ln_limdiahsb, ln_limctl, iiceprt, jiceprt
[5123]	366	!!-------------------------------------------------------------------
[7646]	367	!
[5123]	368	REWIND( numnam_ice_ref ) ! Namelist namicerun in reference namelist : Parameters for ice
	369	READ ( numnam_ice_ref, namicerun, IOSTAT = ios, ERR = 901)
[7646]	370	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in reference namelist', lwp )
	371
[5123]	372	REWIND( numnam_ice_cfg ) ! Namelist namicerun in configuration namelist : Parameters for ice
	373	READ ( numnam_ice_cfg, namicerun, IOSTAT = ios, ERR = 902 )
[7646]	374	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicerun in configuration namelist', lwp )
	375	IF(lwm) WRITE ( numoni, namicerun )
[5123]	376	!
[7646]	377	REWIND( numnam_ice_ref ) ! Namelist namicediag in reference namelist : Parameters for ice
	378	READ ( numnam_ice_ref, namicediag, IOSTAT = ios, ERR = 903)
	379	903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in reference namelist', lwp )
	380
	381	REWIND( numnam_ice_cfg ) ! Namelist namicediag in configuration namelist : Parameters for ice
	382	READ ( numnam_ice_cfg, namicediag, IOSTAT = ios, ERR = 904 )
	383	904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicediag in configuration namelist', lwp )
	384	IF(lwm) WRITE ( numoni, namicediag )
	385	!
[5123]	386	IF(lwp) THEN ! control print
	387	WRITE(numout,*)
[7646]	388	WRITE(numout,*) 'lim_run_init : ice share parameters for dynamics/advection/thermo of sea-ice'
[5123]	389	WRITE(numout,*) ' ~~~~~~'
	390	WRITE(numout,*) ' number of ice categories = ', jpl
	391	WRITE(numout,*) ' number of ice layers = ', nlay_i
	392	WRITE(numout,*) ' number of snow layers = ', nlay_s
[6403]	393	WRITE(numout,*) ' maximum ice concentration for NH = ', rn_amax_n
	394	WRITE(numout,*) ' maximum ice concentration for SH = ', rn_amax_s
[7646]	395	WRITE(numout,*) ' Ice thermodynamics (T) or not (F) ln_limthd = ', ln_limthd
	396	WRITE(numout,*) ' Ice dynamics (T) or not (F) ln_limdyn = ', ln_limdyn
	397	WRITE(numout,*) ' (ln_limdyn=T) Ice dynamics switch nn_limdyn = ', nn_limdyn
	398	WRITE(numout,*) ' 2: total'
	399	WRITE(numout,*) ' 1: advection only (no diffusion, no ridging/rafting)'
	400	WRITE(numout,*) ' 0: advection only (as 1 + prescribed velocity, bypass rheology)'
	401	WRITE(numout,*) ' (ln_limdyn=T) prescribed u-vel (case nn_limdyn=0) = ', rn_uice
	402	WRITE(numout,*) ' (ln_limdyn=T) prescribed v-vel (case nn_limdyn=0) = ', rn_vice
	403	WRITE(numout,*)
	404	WRITE(numout,*) '...and ice diagnostics'
	405	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~'
	406	WRITE(numout,*) ' Diagnose online heat/mass/salt budget ln_limdiachk = ', ln_limdiachk
	407	WRITE(numout,*) ' Output heat/mass/salt budget ln_limdiahsb = ', ln_limdiahsb
	408	WRITE(numout,*) ' control prints in ocean.out for (i,j)=(iiceprt,jiceprt) = ', ln_limctl
	409	WRITE(numout,*) ' i-index for control prints (ln_limctl=true) = ', iiceprt
	410	WRITE(numout,*) ' j-index for control prints (ln_limctl=true) = ', jiceprt
[5123]	411	ENDIF
	412	!
	413	! sea-ice timestep and inverse
[7646]	414	rdt_ice = REAL(nn_fsbc) * rdt
[5123]	415	r1_rdtice = 1._wp / rdt_ice
	416
	417	! inverse of nlay_i and nlay_s
	418	r1_nlay_i = 1._wp / REAL( nlay_i, wp )
	419	r1_nlay_s = 1._wp / REAL( nlay_s, wp )
	420	!
[7646]	421	IF( lwp .AND. ln_bdy .AND. ln_limdiachk ) &
[7761]	422	& CALL ctl_warn('online conservation check activated but it does not work with BDY')
[5167]	423	!
[7646]	424	IF( lwp ) WRITE(numout,*) ' ice timestep rdt_ice = ', rdt_ice
	425	!
	426	END SUBROUTINE lim_run_init
[5123]	427
	428
	429	SUBROUTINE lim_itd_init
	430	!!------------------------------------------------------------------
	431	!! * ROUTINE lim_itd_init *
	432	!!
	433	!! ** Purpose : Initializes the ice thickness distribution
	434	!! ** Method : ...
	435	!! ** input : Namelist namiceitd
	436	!!-------------------------------------------------------------------
	437	INTEGER :: ios ! Local integer output status for namelist read
	438	NAMELIST/namiceitd/ nn_catbnd, rn_himean
	439	!
	440	INTEGER :: jl ! dummy loop index
	441	REAL(wp) :: zc1, zc2, zc3, zx1 ! local scalars
	442	REAL(wp) :: zhmax, znum, zden, zalpha !
	443	!!------------------------------------------------------------------
	444	!
	445	REWIND( numnam_ice_ref ) ! Namelist namiceitd in reference namelist : Parameters for ice
[7646]	446	READ ( numnam_ice_ref, namiceitd, IOSTAT = ios, ERR = 905)
	447	905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in reference namelist', lwp )
	448
[5123]	449	REWIND( numnam_ice_cfg ) ! Namelist namiceitd in configuration namelist : Parameters for ice
[7646]	450	READ ( numnam_ice_cfg, namiceitd, IOSTAT = ios, ERR = 906 )
	451	906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceitd in configuration namelist', lwp )
	452	IF(lwm) WRITE ( numoni, namiceitd )
[5123]	453	!
	454	IF(lwp) THEN ! control print
	455	WRITE(numout,*)
[7646]	456	WRITE(numout,*) 'lim_itd_init : Initialization of ice cat distribution '
	457	WRITE(numout,*) '~~~~~~~~~~~~'
	458	WRITE(numout,*) ' shape of ice categories distribution nn_catbnd = ', nn_catbnd
	459	WRITE(numout,*) ' mean ice thickness in the domain (only active if nn_catbnd=2) rn_himean = ', rn_himean
[5123]	460	ENDIF
[6140]	461	!
[5123]	462	!----------------------------------
[7646]	463	!- Thickness categories boundaries
[5123]	464	!----------------------------------
[6140]	465	!
[5123]	466	hi_max(:) = 0._wp
[6140]	467	!
	468	SELECT CASE ( nn_catbnd ) ! type of ice categories distribution
	469	!
	470	CASE (1) !== tanh function (CICE) ==!
[5123]	471	zc1 = 3._wp / REAL( jpl, wp )
	472	zc2 = 10._wp * zc1
	473	zc3 = 3._wp
	474	DO jl = 1, jpl
	475	zx1 = REAL( jl-1, wp ) / REAL( jpl, wp )
	476	hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) )
	477	END DO
[6140]	478	!
	479	CASE (2) !== h^(-alpha) function ==!
	480	zalpha = 0.05_wp
	481	zhmax = 3._wp * rn_himean
[7646]	482	DO jl = 1, jpl
[5123]	483	znum = jpl * ( zhmax+1 )**zalpha
[6140]	484	zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
[5123]	485	hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
	486	END DO
[6140]	487	!
[5123]	488	END SELECT
[6140]	489	!
	490	DO jl = 1, jpl ! mean thickness by category
[5123]	491	hi_mean(jl) = ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp
	492	END DO
[6140]	493	!
	494	hi_max(jpl) = 99._wp ! set to a big value to ensure that all ice is thinner than hi_max(jpl)
	495	!
[7646]	496	IF(lwp) WRITE(numout,*)
	497	IF(lwp) WRITE(numout,*) ' Thickness category boundaries '
	498	IF(lwp) WRITE(numout,*) ' hi_max ', hi_max(0:jpl)
[5123]	499	!
	500	END SUBROUTINE lim_itd_init
	501
[7646]	502
	503	SUBROUTINE ice_lim_flx( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx )
[4990]	504	!!---------------------------------------------------------------------
[5123]	505	!! * ROUTINE ice_lim_flx *
[7646]	506	!!
[4990]	507	!! ** Purpose : update the ice surface boundary condition by averaging and / or
[7646]	508	!! redistributing fluxes on ice categories
[4990]	509	!!
[7646]	510	!! ** Method : average then redistribute
[4990]	511	!!
[7646]	512	!! ** Action :
[4990]	513	!!---------------------------------------------------------------------
[7646]	514	INTEGER , INTENT(in ) :: k_limflx ! =-1 do nothing; =0 average ;
	515	! ! = 1 average and redistribute ; =2 redistribute
	516	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: ptn_ice ! ice surface temperature
[4990]	517	REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: palb_ice ! ice albedo
	518	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqns_ice ! non solar flux
	519	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pqsr_ice ! net solar flux
	520	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdqn_ice ! non solar flux sensitivity
[5407]	521	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pevap_ice ! sublimation
	522	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pdevap_ice ! sublimation sensitivity
[4990]	523	!
[7753]	524	INTEGER :: jl ! dummy loop index
[4990]	525	!
	526	REAL(wp), POINTER, DIMENSION(:,:) :: zalb_m ! Mean albedo over all categories
	527	REAL(wp), POINTER, DIMENSION(:,:) :: ztem_m ! Mean temperature over all categories
	528	!
	529	REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_m ! Mean solar heat flux over all categories
	530	REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_m ! Mean non solar heat flux over all categories
[5407]	531	REAL(wp), POINTER, DIMENSION(:,:) :: z_evap_m ! Mean sublimation over all categories
[4990]	532	REAL(wp), POINTER, DIMENSION(:,:) :: z_dqn_m ! Mean d(qns)/dT over all categories
[5407]	533	REAL(wp), POINTER, DIMENSION(:,:) :: z_devap_m ! Mean d(evap)/dT over all categories
[4990]	534	!!----------------------------------------------------------------------
[6140]	535	!
[4990]	536	IF( nn_timing == 1 ) CALL timing_start('ice_lim_flx')
	537	!
	538	SELECT CASE( k_limflx ) !== averaged on all ice categories ==!
	539	CASE( 0 , 1 )
[5407]	540	CALL wrk_alloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m)
[4990]	541	!
[5410]	542	z_qns_m (:,:) = fice_ice_ave ( pqns_ice (:,:,:) )
	543	z_qsr_m (:,:) = fice_ice_ave ( pqsr_ice (:,:,:) )
	544	z_dqn_m (:,:) = fice_ice_ave ( pdqn_ice (:,:,:) )
	545	z_evap_m (:,:) = fice_ice_ave ( pevap_ice (:,:,:) )
[5407]	546	z_devap_m(:,:) = fice_ice_ave ( pdevap_ice (:,:,:) )
[4990]	547	DO jl = 1, jpl
[7753]	548	pdqn_ice (:,:,jl) = z_dqn_m(:,:)
	549	pdevap_ice(:,:,jl) = z_devap_m(:,:)
[4990]	550	END DO
	551	!
	552	DO jl = 1, jpl
[7753]	553	pqns_ice (:,:,jl) = z_qns_m(:,:)
	554	pqsr_ice (:,:,jl) = z_qsr_m(:,:)
	555	pevap_ice(:,:,jl) = z_evap_m(:,:)
[4990]	556	END DO
	557	!
[5407]	558	CALL wrk_dealloc( jpi,jpj, z_qsr_m, z_qns_m, z_evap_m, z_dqn_m, z_devap_m)
[4990]	559	END SELECT
[6140]	560	!
[4990]	561	SELECT CASE( k_limflx ) !== redistribution on all ice categories ==!
	562	CASE( 1 , 2 )
	563	CALL wrk_alloc( jpi,jpj, zalb_m, ztem_m )
	564	!
[7646]	565	zalb_m(:,:) = fice_ice_ave ( palb_ice (:,:,:) )
	566	ztem_m(:,:) = fice_ice_ave ( ptn_ice (:,:,:) )
[4990]	567	DO jl = 1, jpl
[7753]	568	pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
	569	pevap_ice(:,:,jl) = pevap_ice(:,:,jl) + pdevap_ice(:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) )
	570	pqsr_ice (:,:,jl) = pqsr_ice (:,:,jl) * ( 1._wp - palb_ice(:,:,jl) ) / ( 1._wp - zalb_m(:,:) )
[4990]	571	END DO
	572	!
	573	CALL wrk_dealloc( jpi,jpj, zalb_m, ztem_m )
	574	END SELECT
	575	!
	576	IF( nn_timing == 1 ) CALL timing_stop('ice_lim_flx')
	577	!
	578	END SUBROUTINE ice_lim_flx
[888]	579
[6140]	580
[5167]	581	SUBROUTINE sbc_lim_bef
[5123]	582	!!----------------------------------------------------------------------
[5167]	583	!! * ROUTINE sbc_lim_bef *
[5123]	584	!!
[7646]	585	!! ** purpose : store ice variables at "before" time step
[5123]	586	!!----------------------------------------------------------------------
[7753]	587	a_i_b (:,:,:) = a_i (:,:,:) ! ice area
	588	e_i_b (:,:,:,:) = e_i (:,:,:,:) ! ice thermal energy
	589	v_i_b (:,:,:) = v_i (:,:,:) ! ice volume
	590	v_s_b (:,:,:) = v_s (:,:,:) ! snow volume
	591	e_s_b (:,:,:,:) = e_s (:,:,:,:) ! snow thermal energy
	592	smv_i_b(:,:,:) = smv_i(:,:,:) ! salt content
	593	oa_i_b (:,:,:) = oa_i (:,:,:) ! areal age content
	594	u_ice_b(:,:) = u_ice(:,:)
	595	v_ice_b(:,:) = v_ice(:,:)
	596	!
	597	at_i_b (:,:) = SUM( a_i_b(:,:,:), dim=3 )
[7646]	598
[5167]	599	END SUBROUTINE sbc_lim_bef
[921]	600
[6140]	601
[5123]	602	SUBROUTINE sbc_lim_diag0
	603	!!----------------------------------------------------------------------
	604	!! * ROUTINE sbc_lim_diag0 *
[888]	605	!!
[5123]	606	!! ** purpose : set ice-ocean and ice-atm. fluxes to zeros at the beggining
	607	!! of the time step
	608	!!----------------------------------------------------------------------
[7753]	609	sfx (:,:) = 0._wp ;
	610	sfx_bri(:,:) = 0._wp ; sfx_lam(:,:) = 0._wp
	611	sfx_sni(:,:) = 0._wp ; sfx_opw(:,:) = 0._wp
	612	sfx_bog(:,:) = 0._wp ; sfx_dyn(:,:) = 0._wp
	613	sfx_bom(:,:) = 0._wp ; sfx_sum(:,:) = 0._wp
	614	sfx_res(:,:) = 0._wp ; sfx_sub(:,:) = 0._wp
	615	!
	616	wfx_snw(:,:) = 0._wp ; wfx_ice(:,:) = 0._wp
	617	wfx_sni(:,:) = 0._wp ; wfx_opw(:,:) = 0._wp
	618	wfx_bog(:,:) = 0._wp ; wfx_dyn(:,:) = 0._wp
	619	wfx_bom(:,:) = 0._wp ; wfx_sum(:,:) = 0._wp
	620	wfx_res(:,:) = 0._wp ; wfx_sub(:,:) = 0._wp
	621	wfx_spr(:,:) = 0._wp ; wfx_lam(:,:) = 0._wp
	622
	623	hfx_thd(:,:) = 0._wp ;
	624	hfx_snw(:,:) = 0._wp ; hfx_opw(:,:) = 0._wp
	625	hfx_bog(:,:) = 0._wp ; hfx_dyn(:,:) = 0._wp
	626	hfx_bom(:,:) = 0._wp ; hfx_sum(:,:) = 0._wp
	627	hfx_res(:,:) = 0._wp ; hfx_sub(:,:) = 0._wp
	628	hfx_spr(:,:) = 0._wp ; hfx_dif(:,:) = 0._wp
	629	hfx_err(:,:) = 0._wp ; hfx_err_rem(:,:) = 0._wp
	630	hfx_err_dif(:,:) = 0._wp
	631	wfx_err_sub(:,:) = 0._wp
	632	!
	633	afx_tot(:,:) = 0._wp ;
	634	afx_dyn(:,:) = 0._wp ; afx_thd(:,:) = 0._wp
	635	!
	636	diag_heat(:,:) = 0._wp ; diag_smvi(:,:) = 0._wp
	637	diag_vice(:,:) = 0._wp ; diag_vsnw(:,:) = 0._wp
[7646]	638
[7753]	639	tau_icebfr(:,:) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
[7646]	640
[5123]	641	END SUBROUTINE sbc_lim_diag0
[5407]	642
[7646]	643
[4990]	644	FUNCTION fice_cell_ave ( ptab )
	645	!!--------------------------------------------------------------------------
	646	!! * Compute average over categories, for grid cell (ice covered and free ocean)
	647	!!--------------------------------------------------------------------------
	648	REAL (wp), DIMENSION (jpi,jpj) :: fice_cell_ave
	649	REAL (wp), DIMENSION (jpi,jpj,jpl), INTENT (in) :: ptab
	650	INTEGER :: jl ! Dummy loop index
[7646]	651
	652	fice_cell_ave (:,:) = 0._wp
[4990]	653	DO jl = 1, jpl
[5123]	654	fice_cell_ave (:,:) = fice_cell_ave (:,:) + a_i (:,:,jl) * ptab (:,:,jl)
[4990]	655	END DO
[7646]	656
[4990]	657	END FUNCTION fice_cell_ave
[7646]	658
	659
[4990]	660	FUNCTION fice_ice_ave ( ptab )
	661	!!--------------------------------------------------------------------------
	662	!! * Compute average over categories, for ice covered part of grid cell
	663	!!--------------------------------------------------------------------------
	664	REAL (kind=wp), DIMENSION (jpi,jpj) :: fice_ice_ave
	665	REAL (kind=wp), DIMENSION (jpi,jpj,jpl), INTENT(in) :: ptab
[888]	666
[4990]	667	fice_ice_ave (:,:) = 0.0_wp
[5167]	668	WHERE ( at_i (:,:) > 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:)
[4990]	669
	670	END FUNCTION fice_ice_ave
	671
[888]	672	#else
	673	!!----------------------------------------------------------------------
	674	!! Default option Dummy module NO LIM 3.0 sea-ice model
	675	!!----------------------------------------------------------------------
	676	CONTAINS
[7777]	677	!
[7646]	678	SUBROUTINE sbc_ice_lim ( kt, ksbc ) ! Dummy routine
[7777]	679	INTEGER, INTENT(in) :: kt, ksbc
[7646]	680	WRITE(,) 'sbc_ice_lim: You should not have seen this print! error?', kt
[888]	681	END SUBROUTINE sbc_ice_lim
[7777]	682	!
[5123]	683	SUBROUTINE sbc_lim_init ! Dummy routine
	684	END SUBROUTINE sbc_lim_init
[888]	685	#endif
	686
	687	!!======================================================================
	688	END MODULE sbcice_lim

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

Download in other formats: