Context Navigation

limtrp.F90 @ 9270

Last change on this file since 9270 was 7993, checked in by frrh, 7 years ago

Merge in missing revisions 6428:2477 inclusive and 6482 from nemo_v3_6_STABLE
branch. In ptic, this includes the fix for restartability of runoff fields in coupled
models. Evolution of coupled models will therefor be affected.

These changes donot affect evolution of the current stand-alone NEMO-CICE GO6
standard configuration.

Work and testing documented in Met Office GMED ticket 320.

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90 @ 9270

Download in other formats: