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 branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 8158

Last change on this file since 8158 was 8152, checked in by vancop, 7 years ago
SIMIP diagnostics, phase 2, commit#3
Property svn:keywords set to `Id`
File size: 28.4 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	!
[6476]	65	INTEGER :: ji, jj, jk, jm , 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	!
[5134]	70	REAL(wp), POINTER, DIMENSION(:,:) :: zsm
	71	REAL(wp), POINTER, DIMENSION(:,:,:) :: z0ice, z0snw, z0ai, z0es , z0smi , z0oi
	72	REAL(wp), POINTER, DIMENSION(:,:,:) :: z0opw
	73	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: z0ei
[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
[6476]	77	REAL(wp), POINTER, DIMENSION(:,:,:) :: zhdfptab
[5123]	78	REAL(wp) :: zdv, zvi, zvs, zsmv, zes, zei
	79	REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
[2715]	80	!!---------------------------------------------------------------------
[6476]	81	INTEGER :: ihdf_vars = 6 !!Number of variables in which we apply horizontal diffusion
	82	!! inside limtrp for each ice category , not counting the
	83	!! variables corresponding to ice_layers
	84	!!---------------------------------------------------------------------
[4161]	85	IF( nn_timing == 1 ) CALL timing_start('limtrp')
[825]	86
[5180]	87	CALL wrk_alloc( jpi,jpj, zsm, zatold, zeiold, zesold )
	88	CALL wrk_alloc( jpi,jpj,jpl, z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
	89	CALL wrk_alloc( jpi,jpj,1, z0opw )
	90	CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei )
	91	CALL wrk_alloc( jpi,jpj,jpl, zhimax, zviold, zvsold, zsmvold )
[6476]	92	CALL wrk_alloc( jpi,jpj,jpl*(ihdf_vars + nlay_i)+1,zhdfptab)
[825]	93
[2715]	94	IF( numit == nstart .AND. lwp ) THEN
	95	WRITE(numout,*)
	96	IF( ln_limdyn ) THEN ; WRITE(numout,*) 'lim_trp : Ice transport '
	97	ELSE ; WRITE(numout,*) 'lim_trp : No ice advection as ln_limdyn = ', ln_limdyn
	98	ENDIF
	99	WRITE(numout,*) '~~~~~~~~~~~~'
[5123]	100	ncfl = 0 ! nb of time step with CFL > 1/2
[2715]	101	ENDIF
[5123]	102
	103	zsm(:,:) = e12t(:,:)
[8150]	104
[2715]	105	! !-------------------------------------!
	106	IF( ln_limdyn ) THEN ! Advection of sea ice properties !
	107	! !-------------------------------------!
[4688]	108
	109	! conservation test
[5123]	110	IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[4688]	111
[5123]	112	! mass and salt flux init
[4161]	113	zviold(:,:,:) = v_i(:,:,:)
[5123]	114	zvsold(:,:,:) = v_s(:,:,:)
	115	zsmvold(:,:,:) = smv_i(:,:,:)
	116	zeiold(:,:) = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )
	117	zesold(:,:) = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )
[825]	118
[8152]	119	! SIMIP diags init
	120	diag_dmtx_dyn(:,:) = 0._wp ; diag_dmty_dyn(:,:) = 0._wp
	121
[5123]	122	!--- Thickness correction init. -------------------------------
	123	zatold(:,:) = SUM( a_i(:,:,:), dim=3 )
[5167]	124	DO jl = 1, jpl
	125	DO jj = 1, jpj
	126	DO ji = 1, jpi
	127	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
	128	ht_i (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
	129	ht_s (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
	130	END DO
	131	END DO
	132	END DO
[4161]	133	!---------------------------------------------------------------------
[5167]	134	! Record max of the surrounding ice thicknesses for correction
[4161]	135	! in case advection creates ice too thick.
	136	!---------------------------------------------------------------------
[5134]	137	zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:)
[4161]	138	DO jl = 1, jpl
	139	DO jj = 2, jpjm1
	140	DO ji = 2, jpim1
[5134]	141	zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) )
[4161]	142	END DO
	143	END DO
	144	CALL lbc_lnk(zhimax(:,:,jl),'T',1.)
	145	END DO
	146
[5123]	147	!=============================!
	148	!== Prather scheme ==!
	149	!=============================!
	150
	151	! If ice drift field is too fast, use an appropriate time step for advection.
	152	zcfl = MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) ) ! CFL test for stability
	153	zcfl = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
	154	IF(lk_mpp ) CALL mpp_max( zcfl )
	155
	156	IF( zcfl > 0.5 ) THEN ; initad = 2 ; zusnit = 0.5_wp
	157	ELSE ; initad = 1 ; zusnit = 1.0_wp
	158	ENDIF
	159
	160	IF( zcfl > 0.5_wp .AND. lwp ) ncfl = ncfl + 1
[5202]	161	!! IF( lwp ) THEN
	162	!! IF( ncfl > 0 ) THEN
	163	!! WRITE(cltmp,'(i6.1)') ncfl
	164	!! CALL ctl_warn( 'lim_trp: ncfl= ', TRIM(cltmp), 'advective ice time-step using a split in sub-time-step ')
	165	!! ELSE
	166	!! ! WRITE(numout,*) 'lim_trp : CFL criterion for ice advection is always smaller than 1/2 '
	167	!! ENDIF
	168	!! ENDIF
[5123]	169
[825]	170	!-------------------------
[2715]	171	! transported fields
[825]	172	!-------------------------
[5134]	173	z0opw(:,:,1) = ato_i(:,:) * e12t(:,:) ! Open water area
[2715]	174	DO jl = 1, jpl
[5134]	175	z0snw (:,:,jl) = v_s (:,:,jl) * e12t(:,:) ! Snow volume
	176	z0ice(:,:,jl) = v_i (:,:,jl) * e12t(:,:) ! Ice volume
	177	z0ai (:,:,jl) = a_i (:,:,jl) * e12t(:,:) ! Ice area
	178	z0smi (:,:,jl) = smv_i(:,:,jl) * e12t(:,:) ! Salt content
	179	z0oi (:,:,jl) = oa_i (:,:,jl) * e12t(:,:) ! Age content
	180	z0es (:,:,jl) = e_s (:,:,1,jl) * e12t(:,:) ! Snow heat content
[6476]	181	DO jk = 1, nlay_i
[5134]	182	z0ei (:,:,jk,jl) = e_i (:,:,jk,jl) * e12t(:,:) ! Ice heat content
[5123]	183	END DO
[825]	184	END DO
	185
[921]	186
[2715]	187	IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN !== odd ice time step: adv_x then adv_y ==!
[5123]	188	DO jt = 1, initad
[5134]	189	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area
[5123]	190	& sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[5134]	191	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:), &
[5123]	192	& sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[825]	193	DO jl = 1, jpl
[8152]	194
[7506]	195	! SIMIP mass transport diags
[8150]	196	diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) - ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) )
[8152]	197
[5134]	198	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume ---
[2715]	199	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[5134]	200	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---
[2715]	201	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[8152]	202
[8150]	203	diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) + ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) )
[7506]	204
[8150]	205	diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) - ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) )
[8152]	206
[7506]	207	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume
	208	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
	209	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume
[2715]	210	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[8152]	211
[8150]	212	diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) + ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) )
[8152]	213
[5134]	214	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---
[2715]	215	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
[5134]	216	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), &
[2715]	217	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
[5134]	218	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---
[2715]	219	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[5134]	220	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), &
[2715]	221	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[5134]	222	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---
[2715]	223	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[5134]	224	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), &
[2715]	225	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[5134]	226	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---
[2715]	227	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[5134]	228	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), &
[2715]	229	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[5123]	230	DO jk = 1, nlay_i !--- ice heat contents ---
[5134]	231	CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), &
[4870]	232	& sxxe(:,:,jk,jl), sye (:,:,jk,jl), &
	233	& syye(:,:,jk,jl), sxye(:,:,jk,jl) )
[5134]	234	CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), &
[4870]	235	& sxxe(:,:,jk,jl), sye (:,:,jk,jl), &
	236	& syye(:,:,jk,jl), sxye(:,:,jk,jl) )
[825]	237	END DO
	238	END DO
	239	END DO
	240	ELSE
[5123]	241	DO jt = 1, initad
[5134]	242	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0opw (:,:,1), sxopw(:,:), & !--- ice open water area
[5123]	243	& sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[5134]	244	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0opw (:,:,1), sxopw(:,:), &
[5123]	245	& sxxopw(:,:) , syopw(:,:), syyopw(:,:), sxyopw(:,:) )
[825]	246	DO jl = 1, jpl
[8152]	247
[8150]	248	! SIMIP mass transport diags
	249	diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) - ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) )
[5134]	250	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume ---
[2715]	251	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[5134]	252	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---
[2715]	253	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[8150]	254	diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) + ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) )
[7506]	255
[8150]	256	diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) - ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) )
[8152]	257	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl), & !--- ice volume ---
[8150]	258	& sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl) )
[8152]	259	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl), & !--- snow volume ---
[2715]	260	& sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl) )
[8150]	261	diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) + ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) )
[7506]	262
[5134]	263	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), & !--- ice salinity ---
[2715]	264	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
[5134]	265	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl), &
[2715]	266	& sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl) )
[5134]	267	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), & !--- ice age ---
[2715]	268	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[5134]	269	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0oi (:,:,jl), sxage(:,:,jl), &
[2715]	270	& sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl) )
[5134]	271	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), & !--- ice concentrations ---
[2715]	272	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[5134]	273	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ai (:,:,jl), sxa (:,:,jl), &
[2715]	274	& sxxa (:,:,jl), sya (:,:,jl), syya (:,:,jl), sxya (:,:,jl) )
[5134]	275	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), & !--- snow heat contents ---
[2715]	276	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[5134]	277	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0es (:,:,jl), sxc0 (:,:,jl), &
[2715]	278	& sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl) )
[4870]	279	DO jk = 1, nlay_i !--- ice heat contents ---
[5134]	280	CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), &
[4870]	281	& sxxe(:,:,jk,jl), sye (:,:,jk,jl), &
	282	& syye(:,:,jk,jl), sxye(:,:,jk,jl) )
[5134]	283	CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ei(:,:,jk,jl), sxe (:,:,jk,jl), &
[4870]	284	& sxxe(:,:,jk,jl), sye (:,:,jk,jl), &
	285	& syye(:,:,jk,jl), sxye(:,:,jk,jl) )
[825]	286	END DO
	287	END DO
	288	END DO
	289	ENDIF
	290
[7506]	291	! SIMIP diags
[8150]	292	diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) / rdt_ice
	293	diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) / rdt_ice
[7506]	294
[825]	295	!-------------------------------------------
	296	! Recover the properties from their contents
	297	!-------------------------------------------
[5134]	298	ato_i(:,:) = z0opw(:,:,1) * r1_e12t(:,:)
[825]	299	DO jl = 1, jpl
[5134]	300	v_i (:,:,jl) = z0ice(:,:,jl) * r1_e12t(:,:)
	301	v_s (:,:,jl) = z0snw(:,:,jl) * r1_e12t(:,:)
	302	smv_i(:,:,jl) = z0smi(:,:,jl) * r1_e12t(:,:)
	303	oa_i (:,:,jl) = z0oi (:,:,jl) * r1_e12t(:,:)
	304	a_i (:,:,jl) = z0ai (:,:,jl) * r1_e12t(:,:)
	305	e_s (:,:,1,jl) = z0es (:,:,jl) * r1_e12t(:,:)
[5123]	306	DO jk = 1, nlay_i
[5134]	307	e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:)
[5123]	308	END DO
[825]	309	END DO
	310
[5123]	311	at_i(:,:) = a_i(:,:,1) ! total ice fraction
	312	DO jl = 2, jpl
	313	at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
	314	END DO
	315
[921]	316	!------------------------------------------------------------------------------!
[5123]	317	! Diffusion of Ice fields
[921]	318	!------------------------------------------------------------------------------!
[7597]	319	IF( nn_ahi0 /= -1 ) THEN
	320
	321	! --- Prepare diffusion for variables with categories --- !
	322	! mask eddy diffusivity coefficient at ocean U- and V-points
	323	jm=1
	324	DO jl = 1, jpl
	325	DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row
	326	DO ji = 1 , fs_jpim1 ! vector opt.
	327	pahu3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji ,jj, jl ) ) ) ) &
	328	& * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj, jl ) ) ) ) * ahiu(ji,jj)
	329	pahv3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji, jj, jl ) ) ) ) &
	330	& * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji, jj+1,jl ) ) ) ) * ahiv(ji,jj)
	331	END DO
[6476]	332	END DO
[7597]	333
	334	zhdfptab(:,:,jm)= a_i (:,:, jl); jm = jm + 1
	335	zhdfptab(:,:,jm)= v_i (:,:, jl); jm = jm + 1
	336	zhdfptab(:,:,jm)= v_s (:,:, jl); jm = jm + 1
	337	zhdfptab(:,:,jm)= smv_i(:,:, jl); jm = jm + 1
	338	zhdfptab(:,:,jm)= oa_i (:,:, jl); jm = jm + 1
	339	zhdfptab(:,:,jm)= e_s (:,:,1,jl); jm = jm + 1
	340	! Sample of adding more variables to apply lim_hdf (ihdf_vars must be increased)
	341	! zhdfptab(:,:,jm) = variable_1 (:,:,1,jl); jm = jm + 1
	342	! zhdfptab(:,:,jm) = variable_2 (:,:,1,jl); jm = jm + 1
	343	DO jk = 1, nlay_i
	344	zhdfptab(:,:,jm)=e_i(:,:,jk,jl); jm= jm+1
	345	END DO
[6476]	346	END DO
[7597]	347	!
	348	! --- Prepare diffusion for open water area --- !
	349	! mask eddy diffusivity coefficient at ocean U- and V-points
	350	DO jj = 1, jpjm1 ! NB: has not to be defined on jpj line and jpi row
	351	DO ji = 1 , fs_jpim1 ! vector opt.
	352	pahu3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji ,jj) ) ) ) &
	353	& * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
	354	pahv3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj ) ) ) ) &
	355	& * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
[825]	356	END DO
	357	END DO
[7597]	358	!
	359	zhdfptab(:,:,jm)= ato_i (:,:);
[825]	360
[7597]	361	! --- Apply diffusion --- !
	362	CALL lim_hdf( zhdfptab, ihdf_vars, jpl, nlay_i)
	363
	364	! --- Recover properties --- !
	365	jm=1
	366	DO jl = 1, jpl
	367	a_i (:,:, jl) = zhdfptab(:,:,jm); jm = jm + 1
	368	v_i (:,:, jl) = zhdfptab(:,:,jm); jm = jm + 1
	369	v_s (:,:, jl) = zhdfptab(:,:,jm); jm = jm + 1
	370	smv_i(:,:, jl) = zhdfptab(:,:,jm); jm = jm + 1
	371	oa_i (:,:, jl) = zhdfptab(:,:,jm); jm = jm + 1
	372	e_s (:,:,1,jl) = zhdfptab(:,:,jm); jm = jm + 1
	373	! Sample of adding more variables to apply lim_hdf---------
	374	! variable_1 (:,:,1,jl) = zhdfptab(:,:, jm ) ; jm + 1
	375	! variable_2 (:,:,1,jl) = zhdfptab(:,:, jm ) ; jm + 1
	376	DO jk = 1, nlay_i
	377	e_i(:,:,jk,jl) = zhdfptab(:,:,jm);jm= jm + 1
[5123]	378	END DO
	379	END DO
[7597]	380	ato_i (:,:) = zhdfptab(:,:,jm)
[4688]	381
[7597]	382	ENDIF
	383
[5167]	384	! --- diags ---
	385	DO jj = 1, jpj
	386	DO ji = 1, jpi
	387	diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice
	388	diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice
	389
	390	diag_trp_vi (ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice
	391	diag_trp_vs (ji,jj) = SUM( v_s(ji,jj,:) - zvsold(ji,jj,:) ) * r1_rdtice
	392	diag_trp_smv(ji,jj) = SUM( smv_i(ji,jj,:) - zsmvold(ji,jj,:) ) * r1_rdtice
	393	END DO
	394	END DO
	395
[5123]	396	! zap small areas
	397	CALL lim_var_zapsmall
	398
	399	!--- Thickness correction in case too high --------------------------------------------------------
[4161]	400	DO jl = 1, jpl
	401	DO jj = 1, jpj
	402	DO ji = 1, jpi
	403
	404	IF ( v_i(ji,jj,jl) > 0._wp ) THEN
[5134]	405
[5167]	406	rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
	407	ht_i (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
	408	ht_s (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
	409
[4688]	410	zvi = v_i (ji,jj,jl)
	411	zvs = v_s (ji,jj,jl)
	412	zsmv = smv_i(ji,jj,jl)
	413	zes = e_s (ji,jj,1,jl)
[4990]	414	zei = SUM( e_i(ji,jj,1:nlay_i,jl) )
[5134]	415
	416	zdv = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl)
	417
	418	IF ( ( zdv > 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. &
[5167]	419	& ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN
[5134]	420
	421	rswitch = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) )
	422	a_i(ji,jj,jl) = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 )
	423
	424	! small correction due to *rswitch for a_i
	425	v_i (ji,jj,jl) = rswitch * v_i (ji,jj,jl)
	426	v_s (ji,jj,jl) = rswitch * v_s (ji,jj,jl)
	427	smv_i(ji,jj,jl) = rswitch * smv_i(ji,jj,jl)
	428	e_s(ji,jj,1,jl) = rswitch * e_s(ji,jj,1,jl)
	429	e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
	430
	431	! Update mass fluxes
	432	wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice
	433	wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
	434	sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice
	435	hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
	436	hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0
	437
[4161]	438	ENDIF
	439
	440	ENDIF
[5123]	441
[825]	442	END DO
	443	END DO
	444	END DO
[4688]	445	! -------------------------------------------------
[5123]	446
	447	!--------------------------------------
	448	! Impose a_i < amax in mono-category
	449	!--------------------------------------
	450	!
	451	IF ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) THEN ! simple conservative piling, comparable with LIM2
	452	DO jj = 1, jpj
	453	DO ji = 1, jpi
[6311]	454	a_i(ji,jj,1) = MIN( a_i(ji,jj,1), rn_amax_2d(ji,jj) )
[5123]	455	END DO
	456	END DO
	457	ENDIF
[825]	458
[4990]	459	! --- agglomerate variables -----------------
[4688]	460	vt_i (:,:) = 0._wp
	461	vt_s (:,:) = 0._wp
	462	at_i (:,:) = 0._wp
[825]	463	DO jl = 1, jpl
	464	DO jj = 1, jpj
	465	DO ji = 1, jpi
[5134]	466	vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl)
	467	vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl)
	468	at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl)
[4688]	469	END DO
	470	END DO
	471	END DO
[825]	472
[5134]	473	! --- open water = 1 if at_i=0 --------------------------------
[4161]	474	DO jj = 1, jpj
	475	DO ji = 1, jpi
[4990]	476	rswitch = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
[5123]	477	ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
[4161]	478	END DO
[4688]	479	END DO
[4161]	480
[4688]	481	! conservation test
	482	IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
[4161]	483
[825]	484	ENDIF
	485
[5123]	486	! -------------------------------------------------
	487	! control prints
	488	! -------------------------------------------------
[5128]	489	IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )
[2715]	490	!
[5180]	491	CALL wrk_dealloc( jpi,jpj, zsm, zatold, zeiold, zesold )
	492	CALL wrk_dealloc( jpi,jpj,jpl, z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
	493	CALL wrk_dealloc( jpi,jpj,1, z0opw )
	494	CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei )
	495	CALL wrk_dealloc( jpi,jpj,jpl, zviold, zvsold, zhimax, zsmvold )
[6476]	496	CALL wrk_dealloc( jpi,jpj,jpl*(ihdf_vars+nlay_i)+1,zhdfptab)
[5123]	497	!
[4161]	498	IF( nn_timing == 1 ) CALL timing_stop('limtrp')
[5123]	499
[825]	500	END SUBROUTINE lim_trp
	501
	502	#else
	503	!!----------------------------------------------------------------------
	504	!! Default option Empty Module No sea-ice model
	505	!!----------------------------------------------------------------------
	506	CONTAINS
	507	SUBROUTINE lim_trp ! Empty routine
	508	END SUBROUTINE lim_trp
	509	#endif
	510	!!======================================================================
	511	END MODULE limtrp
[6476]	512

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

Download in other formats: