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.

limtrp.F90 in trunk/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 5128

Last change on this file since 5128 was 5128, checked in by clem, 9 years ago
LIM3: put debug (i,j) points in the namelist and change their names
Property svn:keywords set to `Id`
File size: 26.1 KB

Rev	Line
[825]	1	MODULE limtrp
	2	!!======================================================================
	3	!! * MODULE limtrp *
	4	!! LIM transport ice model : sea-ice advection/diffusion
	5	!!======================================================================
[2715]	6	!! History : LIM-2 ! 2000-01 (M.A. Morales Maqueda, H. Goosse, and T. Fichefet) Original code
	7	!! 3.0 ! 2005-11 (M. Vancoppenolle) Multi-layer sea ice, salinity variations
	8	!! 4.0 ! 2011-02 (G. Madec) dynamical allocation
	9	!!----------------------------------------------------------------------
[825]	10	#if defined key_lim3
	11	!!----------------------------------------------------------------------
[834]	12	!! 'key_lim3' LIM3 sea-ice model
[825]	13	!!----------------------------------------------------------------------
	14	!! lim_trp : advection/diffusion process of sea ice
	15	!!----------------------------------------------------------------------
[3625]	16	USE phycst ! physical constant
	17	USE dom_oce ! ocean domain
	18	USE sbc_oce ! ocean surface boundary condition
	19	USE dom_ice ! ice domain
	20	USE ice ! ice variables
	21	USE limadv ! ice advection
	22	USE limhdf ! ice horizontal diffusion
[5123]	23	USE limvar !
	24	!
[3625]	25	USE in_out_manager ! I/O manager
	26	USE lbclnk ! lateral boundary conditions -- MPP exchanges
	27	USE lib_mpp ! MPP library
	28	USE wrk_nemo ! work arrays
	29	USE prtctl ! Print control
	30	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[5123]	31	USE timing ! Timing
[4688]	32	USE limcons ! conservation tests
[5123]	33	USE limctl ! control prints
[825]	34
	35	IMPLICIT NONE
	36	PRIVATE
	37
[5123]	38	PUBLIC lim_trp ! called by sbcice_lim
[825]	39
[5123]	40	INTEGER :: ncfl ! number of ice time step with CFL>1/2
	41
[825]	42	!! * Substitution
	43	# include "vectopt_loop_substitute.h90"
	44	!!----------------------------------------------------------------------
[4161]	45	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
[1156]	46	!! $Id$
[2715]	47	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[825]	48	!!----------------------------------------------------------------------
	49	CONTAINS
	50
[921]	51	SUBROUTINE lim_trp( kt )
[825]	52	!!-------------------------------------------------------------------
	53	!! * ROUTINE lim_trp *
	54	!!
	55	!! ** purpose : advection/diffusion process of sea ice
	56	!!
	57	!! ** method : variables included in the process are scalar,
	58	!! other values are considered as second order.
	59	!! For advection, a second order Prather scheme is used.
	60	!!
	61	!! ** action :
	62	!!---------------------------------------------------------------------
[5123]	63	INTEGER, INTENT(in) :: kt ! number of iteration
[2715]	64	!
[5123]	65	INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices
[2715]	66	INTEGER :: initad ! number of sub-timestep for the advection
[4990]	67	REAL(wp) :: zcfl , zusnit ! - -
[5123]	68	CHARACTER(len=80) :: cltmp
[2715]	69	!
[5123]	70	REAL(wp), POINTER, DIMENSION(:,:) :: zsm, zs0at
	71	REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e
	72	REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ow
[3294]	73	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zs0e
[5123]	74	REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold, zsmvold ! old ice volume...
	75	REAL(wp), POINTER, DIMENSION(:,:,:) :: zhimax ! old ice thickness
	76	REAL(wp), POINTER, DIMENSION(:,:) :: zatold, zeiold, zesold ! old concentration, enthalpies
	77	REAL(wp) :: zdv, zvi, zvs, zsmv, zes, zei
	78	REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
[2715]	79	!!---------------------------------------------------------------------
[4161]	80	IF( nn_timing == 1 ) CALL timing_start('limtrp')
[825]	81
[5123]	82	CALL wrk_alloc( jpi,jpj, zsm, zs0at, zatold, zeiold, zesold )
	83	CALL wrk_alloc( jpi,jpj,jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e )
	84	CALL wrk_alloc( jpi,jpj,1, zs0ow )
	85	CALL wrk_alloc( jpi,jpj,nlay_i+1,jpl, zs0e )
	86	CALL wrk_alloc( jpi,jpj,jpl, zhimax, zviold, zvsold, zsmvold )
[825]	87
[2715]	88	IF( numit == nstart .AND. lwp ) THEN
	89	WRITE(numout,*)
	90	IF( ln_limdyn ) THEN ; WRITE(numout,*) 'lim_trp : Ice transport '
	91	ELSE ; WRITE(numout,*) 'lim_trp : No ice advection as ln_limdyn = ', ln_limdyn
	92	ENDIF
	93	WRITE(numout,*) '~~~~~~~~~~~~'
[5123]	94	ncfl = 0 ! nb of time step with CFL > 1/2
[2715]	95	ENDIF
[5123]	96
	97	zsm(:,:) = e12t(:,:)
[2715]	98
	99	! !-------------------------------------!
	100	IF( ln_limdyn ) THEN ! Advection of sea ice properties !
	101	! !-------------------------------------!
[4688]	102
	103	! conservation test
[5123]	104	IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[4688]	105
[5123]	106	! mass and salt flux init
[4161]	107	zviold(:,:,:) = v_i(:,:,:)
[5123]	108	zvsold(:,:,:) = v_s(:,:,:)
	109	zsmvold(:,:,:) = smv_i(:,:,:)
	110	zeiold(:,:) = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )
	111	zesold(:,:) = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )
[825]	112
[5123]	113	!--- Thickness correction init. -------------------------------
[4161]	114	CALL lim_var_glo2eqv
[5123]	115	zatold(:,:) = SUM( a_i(:,:,:), dim=3 )
[4161]	116	!---------------------------------------------------------------------
	117	! Record max of the surrounding ice thicknesses for correction in limupdate
	118	! in case advection creates ice too thick.
	119	!---------------------------------------------------------------------
	120	zhimax(:,:,:) = ht_i(:,:,:)
	121	DO jl = 1, jpl
	122	DO jj = 2, jpjm1
	123	DO ji = 2, jpim1
	124	zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) )
	125	END DO
	126	END DO
	127	CALL lbc_lnk(zhimax(:,:,jl),'T',1.)
	128	END DO
	129
[5123]	130	!=============================!
	131	!== Prather scheme ==!
	132	!=============================!
	133
	134	! If ice drift field is too fast, use an appropriate time step for advection.
	135	zcfl = MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) ! CFL test for stability
	136	zcfl = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
	137	IF(lk_mpp ) CALL mpp_max( zcfl )
	138
	139	IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp
	140	ELSE ; initad = 1 ; zusnit = 1.0_wp
	141	ENDIF
	142
	143	IF( zcfl > 0.5_wp .AND. lwp ) ncfl = ncfl + 1
	144	IF( numit == nlast .AND. lwp ) THEN
	145	IF( ncfl > 0 ) THEN
	146	WRITE(cltmp,'(i6.1)') ncfl
	147	CALL ctl_stop('STOP',TRIM(cltmp) )
	148	CALL ctl_warn( 'lim_trp: ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ')
	149	ELSE
	150	WRITE(numout,*) 'lim_trp : CFL criteria for ice advection is always smaller than 1/2 '
	151	ENDIF
	152	ENDIF
	153
[825]	154	!-------------------------
[2715]	155	! transported fields
[825]	156	!-------------------------
[5123]	157	zs0ow(:,:,1) = ato_i(:,:) * e12t(:,:) ! Open water area
[2715]	158	DO jl = 1, jpl
[5123]	159	zs0sn (:,:,jl) = v_s (:,:,jl) * e12t(:,:) ! Snow volume
	160	zs0ice(:,:,jl) = v_i (:,:,jl) * e12t(:,:) ! Ice volume
	161	zs0a (:,:,jl) = a_i (:,:,jl) * e12t(:,:) ! Ice area
	162	zs0sm (:,:,jl) = smv_i(:,:,jl) * e12t(:,:) ! Salt content
	163	zs0oi (:,:,jl) = oa_i (:,:,jl) * e12t(:,:) ! Age content
	164	zs0c0 (:,:,jl) = e_s (:,:,1,jl) * e12t(:,:) ! Snow heat content
	165	DO jk = 1, nlay_i
	166	zs0e (:,:,jk,jl) = e_i (:,:,jk,jl) * e12t(:,:) ! Ice heat content
	167	END DO
[825]	168	END DO
	169
[921]	170
[2715]	171	IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==!
[5123]	172	DO jt = 1, initad
	173	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0ow (:,:,1), sxopw(:,:), & !--- ice open water area
	174	& sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )
	175	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ow (:,:,1), sxopw(:,:), &
	176	& sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[825]	177	DO jl = 1, jpl
[5123]	178	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---
[2715]	179	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[4688]	180	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl), &
[2715]	181	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[5123]	182	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---
[2715]	183	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[4688]	184	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), &
[2715]	185	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[5123]	186	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---
[2715]	187	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
[4688]	188	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), &
[2715]	189	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
[5123]	190	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---
[2715]	191	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[4688]	192	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl), &
[2715]	193	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[5123]	194	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---
[2715]	195	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[4688]	196	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0a (:,:,jl), sxa (:,:,jl), &
[2715]	197	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[5123]	198	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---
[2715]	199	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[4688]	200	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), &
[2715]	201	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[5123]	202	DO jk = 1, nlay_i !--- ice heat contents ---
	203	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), &
[4870]	204	& sxxe(:,:,jk,jl), sye (:,:,jk,jl), &
	205	& syye(:,:,jk,jl), sxye(:,:,jk,jl) )
	206	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), &
	207	& sxxe(:,:,jk,jl), sye (:,:,jk,jl), &
	208	& syye(:,:,jk,jl), sxye(:,:,jk,jl) )
[825]	209	END DO
	210	END DO
	211	END DO
	212	ELSE
[5123]	213	DO jt = 1, initad
	214	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0ow (:,:,1), sxopw(:,:), & !--- ice open water area
	215	& sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )
	216	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ow (:,:,1), sxopw(:,:), &
	217	& sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[825]	218	DO jl = 1, jpl
[5123]	219	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl), & !--- ice volume ---
[2715]	220	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[4688]	221	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0ice(:,:,jl), sxice(:,:,jl), &
[2715]	222	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[5123]	223	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---
[2715]	224	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[4688]	225	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sn (:,:,jl), sxsn (:,:,jl), &
[2715]	226	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[5123]	227	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---
[2715]	228	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
[4688]	229	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0sm (:,:,jl), sxsal(:,:,jl), &
[2715]	230	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
	231
[5123]	232	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---
[2715]	233	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[4688]	234	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0oi (:,:,jl), sxage(:,:,jl), &
[2715]	235	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[5123]	236	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0a (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---
[2715]	237	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[4688]	238	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0a (:,:,jl), sxa (:,:,jl), &
[2715]	239	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[5123]	240	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---
[2715]	241	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[4688]	242	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl), &
[2715]	243	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[4870]	244	DO jk = 1, nlay_i !--- ice heat contents ---
[5123]	245	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), &
[4870]	246	& sxxe(:,:,jk,jl), sye (:,:,jk,jl), &
	247	& syye(:,:,jk,jl), sxye(:,:,jk,jl) )
	248	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl), &
	249	& sxxe(:,:,jk,jl), sye (:,:,jk,jl), &
	250	& syye(:,:,jk,jl), sxye(:,:,jk,jl) )
[825]	251	END DO
	252	END DO
	253	END DO
	254	ENDIF
	255
	256	!-------------------------------------------
	257	! Recover the properties from their contents
	258	!-------------------------------------------
[5123]	259	ato_i(:,:) = zs0ow(:,:,1) * r1_e12t(:,:)
[825]	260	DO jl = 1, jpl
[5123]	261	v_i (:,:,jl) = zs0ice(:,:,jl) * r1_e12t(:,:)
	262	v_s (:,:,jl) = zs0sn (:,:,jl) * r1_e12t(:,:)
	263	smv_i(:,:,jl) = zs0sm (:,:,jl) * r1_e12t(:,:)
	264	oa_i (:,:,jl) = zs0oi (:,:,jl) * r1_e12t(:,:)
	265	a_i (:,:,jl) = zs0a (:,:,jl) * r1_e12t(:,:)
	266	e_s (:,:,1,jl) = zs0c0 (:,:,jl) * r1_e12t(:,:)
	267	DO jk = 1, nlay_i
	268	e_i (:,:,jk,jl) = zs0e (:,:,jk,jl) * r1_e12t(:,:)
	269	END DO
[825]	270	END DO
	271
[5123]	272	at_i(:,:) = a_i(:,:,1) ! total ice fraction
	273	DO jl = 2, jpl
	274	at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
	275	END DO
	276
[921]	277	!------------------------------------------------------------------------------!
[5123]	278	! Diffusion of Ice fields
[921]	279	!------------------------------------------------------------------------------!
[825]	280
[5123]	281	!
[2715]	282	!--------------------------------
	283	! diffusion of open water area
	284	!--------------------------------
	285	! ! Masked eddy diffusivity coefficient at ocean U- and V-points
	286	DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row
	287	DO ji = 1 , fs_jpim1 ! vector opt.
[5123]	288	pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji ,jj) ) ) ) &
	289	& * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
	290	pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj ) ) ) ) &
	291	& * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
[2715]	292	END DO
	293	END DO
	294	!
[5123]	295	CALL lim_hdf( ato_i (:,:) ) ! Diffusion
[2715]	296
[921]	297	!------------------------------------
[2715]	298	! Diffusion of other ice variables
[921]	299	!------------------------------------
[825]	300	DO jl = 1, jpl
[2715]	301	! ! Masked eddy diffusivity coefficient at ocean U- and V-points
	302	DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row
	303	DO ji = 1 , fs_jpim1 ! vector opt.
[5123]	304	pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji ,jj,jl) ) ) ) &
	305	& * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
	306	pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,jj ,jl) ) ) ) &
	307	& * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
[825]	308	END DO
[921]	309	END DO
[825]	310
[5123]	311	CALL lim_hdf( v_i (:,:, jl) )
	312	CALL lim_hdf( v_s (:,:, jl) )
	313	CALL lim_hdf( smv_i(:,:, jl) )
	314	CALL lim_hdf( oa_i (:,:, jl) )
	315	CALL lim_hdf( a_i (:,:, jl) )
	316	CALL lim_hdf( e_s (:,:,1,jl) )
[825]	317	DO jk = 1, nlay_i
[5123]	318	CALL lim_hdf( e_i(:,:,jk,jl) )
[2715]	319	END DO
	320	END DO
[825]	321
[921]	322	!------------------------------------------------------------------------------!
[5123]	323	! limit ice properties after transport
[921]	324	!------------------------------------------------------------------------------!
[5123]	325	!!gm & cr : MAX should not be active if adv scheme is positive !
[825]	326	DO jl = 1, jpl
	327	DO jj = 1, jpj
	328	DO ji = 1, jpi
[5123]	329	v_s (ji,jj,jl) = MAX( 0._wp, v_s (ji,jj,jl) )
	330	v_i (ji,jj,jl) = MAX( 0._wp, v_i (ji,jj,jl) )
	331	smv_i(ji,jj,jl) = MAX( 0._wp, smv_i(ji,jj,jl) )
	332	oa_i (ji,jj,jl) = MAX( 0._wp, oa_i (ji,jj,jl) )
	333	a_i (ji,jj,jl) = MAX( 0._wp, a_i (ji,jj,jl) )
	334	e_s (ji,jj,1,jl) = MAX( 0._wp, e_s (ji,jj,1,jl) )
[825]	335	END DO
	336	END DO
	337
[4688]	338	DO jk = 1, nlay_i
	339	DO jj = 1, jpj
	340	DO ji = 1, jpi
[5123]	341	e_i(ji,jj,jk,jl) = MAX( 0._wp, e_i(ji,jj,jk,jl) )
	342	END DO
	343	END DO
	344	END DO
	345	END DO
	346	!!gm & cr
[4688]	347
[5123]	348	! zap small areas
	349	CALL lim_var_zapsmall
	350
	351	!--- Thickness correction in case too high --------------------------------------------------------
[4161]	352	CALL lim_var_glo2eqv
	353	DO jl = 1, jpl
	354	DO jj = 1, jpj
	355	DO ji = 1, jpi
	356
	357	IF ( v_i(ji,jj,jl) > 0._wp ) THEN
[4688]	358	zvi = v_i (ji,jj,jl)
	359	zvs = v_s (ji,jj,jl)
	360	zsmv = smv_i(ji,jj,jl)
	361	zes = e_s (ji,jj,1,jl)
[4990]	362	zei = SUM( e_i(ji,jj,1:nlay_i,jl) )
[4688]	363	zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl)
[4161]	364
[4990]	365	rswitch = 1._wp
[5123]	366	IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. &
[4161]	367	& ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN
	368	ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) )
[4990]	369	rswitch = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) )
	370	a_i(ji,jj,jl) = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 )
[4161]	371	ELSE
	372	ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) )
[4990]	373	rswitch = MAX( 0._wp, SIGN( 1._wp, ht_i(ji,jj,jl) - epsi20 ) )
	374	a_i(ji,jj,jl) = rswitch * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi20 )
[4161]	375	ENDIF
	376
[4990]	377	! small correction due to *rswitch for a_i
	378	v_i (ji,jj,jl) = rswitch * v_i (ji,jj,jl)
	379	v_s (ji,jj,jl) = rswitch * v_s (ji,jj,jl)
	380	smv_i(ji,jj,jl) = rswitch * smv_i(ji,jj,jl)
	381	e_s(ji,jj,1,jl) = rswitch * e_s(ji,jj,1,jl)
	382	e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
[4161]	383
	384	! Update mass fluxes
[4688]	385	wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
	386	wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
	387	sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice
[5123]	388	hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
	389	hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0
[4161]	390	ENDIF
[5123]	391
[825]	392	END DO
	393	END DO
	394	END DO
[4688]	395	! -------------------------------------------------
[5123]	396
	397	!--------------------------------------
	398	! Impose a_i < amax in mono-category
	399	!--------------------------------------
	400	!
	401	IF ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) THEN ! simple conservative piling, comparable with LIM2
	402	DO jj = 1, jpj
	403	DO ji = 1, jpi
	404	a_i(ji,jj,1) = MIN( a_i(ji,jj,1), rn_amax )
	405	END DO
	406	END DO
	407	ENDIF
[825]	408
[4688]	409	! --- diags ---
[825]	410	DO jj = 1, jpj
	411	DO ji = 1, jpi
[5123]	412	diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice
	413	diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice
[4990]	414
[5123]	415	diag_trp_vi (ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice
	416	diag_trp_vs (ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice
	417	diag_trp_smv(ji,jj) = SUM( smv_i(ji,jj,:) - zsmvold(ji,jj,:) ) * r1_rdtice
[825]	418	END DO
	419	END DO
	420
[4990]	421	! --- agglomerate variables -----------------
[4688]	422	vt_i (:,:) = 0._wp
	423	vt_s (:,:) = 0._wp
	424	at_i (:,:) = 0._wp
	425	!
[825]	426	DO jl = 1, jpl
	427	DO jj = 1, jpj
	428	DO ji = 1, jpi
[4688]	429	!
	430	vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume
	431	vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume
	432	at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration
	433	END DO
	434	END DO
	435	END DO
	436	! -------------------------------------------------
[825]	437
[4688]	438	! open water
[4161]	439	DO jj = 1, jpj
	440	DO ji = 1, jpi
[4688]	441	! open water = 1 if at_i=0
[4990]	442	rswitch = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
[5123]	443	ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
[4161]	444	END DO
[4688]	445	END DO
[4161]	446
[4688]	447	! conservation test
	448	IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[4161]	449
[825]	450	ENDIF
	451
[5123]	452	CALL lim_var_glo2eqv ! equivalent variables, requested for rafting
[4161]	453
[5123]	454	! -------------------------------------------------
	455	! control prints
	456	! -------------------------------------------------
[5128]	457	IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )
[2715]	458	!
[5123]	459	CALL wrk_dealloc( jpi,jpj, zsm, zs0at, zatold, zeiold, zesold )
	460	CALL wrk_dealloc( jpi,jpj,jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi, zzs0e )
	461	CALL wrk_dealloc( jpi,jpj,1, zs0ow )
	462	CALL wrk_dealloc( jpi,jpj,nlay_i+1,jpl, zs0e )
	463	CALL wrk_dealloc( jpi,jpj,jpl, zviold, zvsold, zhimax, zsmvold )
	464	!
[4161]	465	IF( nn_timing == 1 ) CALL timing_stop('limtrp')
[5123]	466
[825]	467	END SUBROUTINE lim_trp
	468
	469	#else
	470	!!----------------------------------------------------------------------
	471	!! Default option Empty Module No sea-ice model
	472	!!----------------------------------------------------------------------
	473	CONTAINS
	474	SUBROUTINE lim_trp ! Empty routine
	475	END SUBROUTINE lim_trp
	476	#endif
	477	!!======================================================================
	478	END MODULE limtrp

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

Download in other formats: