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.

iscplrst.F90 in branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90 @ 9076

Last change on this file since 9076 was 9019, checked in by timgraham, 7 years ago
Merge of dev_CNRS_2017 into branch
File size: 19.3 KB

Rev	Line
[5790]	1	MODULE iscplrst
	2	!!======================================================================
[7646]	3	!! * MODULE iscplrst *
[5790]	4	!! Ocean forcing: update the restart file in case of ice sheet/ocean coupling
	5	!!=====================================================================
	6	!! History : NEMO ! 2015-01 P. Mathiot: original
	7	!!----------------------------------------------------------------------
	8
	9	!!----------------------------------------------------------------------
[5835]	10	!! iscpl_stp : step management
[5790]	11	!! iscpl_rst_interpol : restart interpolation in case of coupling with ice sheet
	12	!!----------------------------------------------------------------------
[9019]	13	USE oce ! global tra/dyn variable
[5790]	14	USE dom_oce ! ocean space and time domain
	15	USE domwri ! ocean space and time domain
[9019]	16	USE domvvl , ONLY : dom_vvl_interpol
[5790]	17	USE phycst ! physical constants
	18	USE sbc_oce ! surface boundary condition variables
[9019]	19	USE iscplini ! ice sheet coupling: initialisation
	20	USE iscplhsb ! ice sheet coupling: conservation
	21	!
[5790]	22	USE in_out_manager ! I/O manager
	23	USE iom ! I/O module
	24	USE lib_mpp ! MPP library
	25	USE lib_fortran ! MPP library
[5823]	26	USE lbclnk ! communication
[5790]	27
	28	IMPLICIT NONE
	29	PRIVATE
	30
[5823]	31	PUBLIC iscpl_stp ! step management
[5790]	32	!!
	33	!! * Substitutions
[5920]	34	# include "vectopt_loop_substitute.h90"
[5790]	35	!!----------------------------------------------------------------------
	36	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
	37	!! $Id: sbcrnf.F90 4666 2014-06-11 12:52:23Z mathiot $
	38	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
	39	!!----------------------------------------------------------------------
	40	CONTAINS
	41
	42	SUBROUTINE iscpl_stp
	43	!!----------------------------------------------------------------------
	44	!! * ROUTINE iscpl_stp *
	45	!!
	46	!! ** Purpose : compute initialisation
	47	!! compute extrapolation of restart variable un, vn, tsn, sshn (wetting/drying)
	48	!! compute correction term if needed
	49	!!
	50	!!----------------------------------------------------------------------
[5945]	51	INTEGER :: inum0
[9019]	52	REAL(wp), DIMENSION(jpi,jpj) :: zsmask_b
	53	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b, zumask_b, zvmask_b
	54	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t_b , ze3u_b , ze3v_b
	55	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdepw_b
[5945]	56	CHARACTER(20) :: cfile
[5823]	57	!!----------------------------------------------------------------------
[9019]	58	!
	59	! ! get restart variable
[5790]	60	CALL iom_get( numror, jpdom_autoglo, 'tmask' , ztmask_b ) ! need to extrapolate T/S
	61	CALL iom_get( numror, jpdom_autoglo, 'umask' , zumask_b ) ! need to correct barotropic velocity
	62	CALL iom_get( numror, jpdom_autoglo, 'vmask' , zvmask_b ) ! need to correct barotropic velocity
	63	CALL iom_get( numror, jpdom_autoglo, 'smask' , zsmask_b ) ! need to correct barotropic velocity
[6074]	64	CALL iom_get( numror, jpdom_autoglo, 'e3t_n' , ze3t_b(:,:,:) ) ! need to compute temperature correction
	65	CALL iom_get( numror, jpdom_autoglo, 'e3u_n' , ze3u_b(:,:,:) ) ! need to correct barotropic velocity
	66	CALL iom_get( numror, jpdom_autoglo, 'e3v_n' , ze3v_b(:,:,:) ) ! need to correct barotropic velocity
	67	CALL iom_get( numror, jpdom_autoglo, 'gdepw_n', zdepw_b(:,:,:) ) ! need to interpol vertical profile (vvl)
[9019]	68	!
	69	CALL iscpl_init() ! read namelist
	70	! ! Extrapolation/interpolation of modify cell and new cells ... (maybe do it later after domvvl)
[5790]	71	CALL iscpl_rst_interpol( ztmask_b, zumask_b, zvmask_b, zsmask_b, ze3t_b, ze3u_b, ze3v_b, zdepw_b )
[9019]	72	!
	73	IF ( ln_hsb ) THEN ! compute correction if conservation needed
[5790]	74	IF( iscpl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'rst_iscpl : unable to allocate rst_iscpl arrays' )
	75	CALL iscpl_cons(ztmask_b, zsmask_b, ze3t_b, htsc_iscpl, hdiv_iscpl, rdt_iscpl)
	76	END IF
	77
[9019]	78	! ! create a domain file
	79	IF( nn_msh /= 0 .AND. ln_iscpl ) CALL dom_wri
	80	!
[5790]	81	IF ( ln_hsb ) THEN
	82	cfile='correction'
	83	cfile = TRIM( cfile )
	84	CALL iom_open ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib )
	85	CALL iom_rstput( 0, 0, inum0, 'vol_cor', hdiv_iscpl(:,:,:) )
	86	CALL iom_rstput( 0, 0, inum0, 'tem_cor', htsc_iscpl(:,:,:,jp_tem) )
	87	CALL iom_rstput( 0, 0, inum0, 'sal_cor', htsc_iscpl(:,:,:,jp_sal) )
	88	CALL iom_close ( inum0 )
	89	END IF
[9019]	90	!
	91	neuler = 0 ! next step is an euler time step
	92	!
	93	! ! set _b and _n variables equal
[5802]	94	tsb (:,:,:,:) = tsn (:,:,:,:)
[7646]	95	ub (:,:,:) = un (:,:,:)
	96	vb (:,:,:) = vn (:,:,:)
	97	sshb(:,:) = sshn(:,:)
[9019]	98	!
	99	! ! set _b and _n vertical scale factor equal
[6074]	100	e3t_b (:,:,:) = e3t_n (:,:,:)
	101	e3u_b (:,:,:) = e3u_n (:,:,:)
	102	e3v_b (:,:,:) = e3v_n (:,:,:)
[9019]	103	!
[7646]	104	e3uw_b (:,:,:) = e3uw_n (:,:,:)
	105	e3vw_b (:,:,:) = e3vw_n (:,:,:)
	106	gdept_b(:,:,:) = gdept_n(:,:,:)
[6080]	107	gdepw_b(:,:,:) = gdepw_n(:,:,:)
[7646]	108	hu_b (:,:) = hu_n (:,:)
	109	hv_b (:,:) = hv_n (:,:)
	110	r1_hu_b(:,:) = r1_hu_n(:,:)
	111	r1_hv_b(:,:) = r1_hv_n(:,:)
[5790]	112	!
	113	END SUBROUTINE iscpl_stp
[7646]	114
	115
[5790]	116	SUBROUTINE iscpl_rst_interpol (ptmask_b, pumask_b, pvmask_b, psmask_b, pe3t_b, pe3u_b, pe3v_b, pdepw_b)
	117	!!----------------------------------------------------------------------
	118	!! * ROUTINE iscpl_rst_interpol *
	119	!!
	120	!! ** Purpose : compute new tn, sn, un, vn and sshn in case of evolving geometry of ice shelves
	121	!! compute 2d fields of heat, salt and volume correction
	122	!!
	123	!! ** Method : tn, sn : extrapolation from neigbourg cells
	124	!! un, vn : fill with 0 velocity and keep barotropic transport by modifing surface velocity or adjacent velocity
	125	!!----------------------------------------------------------------------
	126	REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: ptmask_b, pumask_b, pvmask_b !! mask before
	127	REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: pe3t_b , pe3u_b , pe3v_b !! scale factor before
[5823]	128	REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: pdepw_b !! depth w before
[5790]	129	REAL(wp), DIMENSION(:,: ), INTENT(in ) :: psmask_b !! mask before
	130	!!
	131	INTEGER :: ji, jj, jk, iz !! loop index
	132	INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
	133	!!
	134	REAL(wp):: summsk, zsum, zsum1, zarea, zsumn, zsumb
	135	REAL(wp):: zdz, zdzm1, zdzp1
	136	!!
[9019]	137	REAL(wp), DIMENSION(jpi,jpj) :: zdmask , zsmask0, zucorr, zbub, zbun, zssh0, zhu1, zde3t
	138	REAL(wp), DIMENSION(jpi,jpj) :: zdsmask, zsmask1, zvcorr, zbvb, zbvn, zssh1, zhv1
	139	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask0, zwmaskn, ztrp
	140	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask1, zwmaskb, ztmp3d
	141	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0
[5945]	142	!!----------------------------------------------------------------------
[9019]	143	!
	144	! ! mask value to be sure
[5790]	145	tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * ptmask_b(:,:,:)
	146	tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * ptmask_b(:,:,:)
[9019]	147	!
	148	! ! compute wmask
[5790]	149	zwmaskn(:,:,1) = tmask (:,:,1)
	150	zwmaskb(:,:,1) = ptmask_b(:,:,1)
	151	DO jk = 2,jpk
	152	zwmaskn(:,:,jk) = tmask (:,:,jk) * tmask (:,:,jk-1)
	153	zwmaskb(:,:,jk) = ptmask_b(:,:,jk) * ptmask_b(:,:,jk-1)
	154	END DO
[9019]	155	!
	156	! ! compute new ssh if we open a full water column (average of the closest neigbourgs)
[5790]	157	sshb (:,:)=sshn(:,:)
	158	zssh0(:,:)=sshn(:,:)
[5802]	159	zsmask0(:,:) = psmask_b(:,:)
	160	zsmask1(:,:) = psmask_b(:,:)
[9019]	161	DO iz = 1, 10 ! need to be tuned (configuration dependent) (OK for ISOMIP+)
[5790]	162	zdsmask(:,:) = ssmask(:,:)-zsmask0(:,:)
[5920]	163	DO jj = 2,jpj-1
	164	DO ji = fs_2, fs_jpim1 ! vector opt.
[5790]	165	jip1=ji+1; jim1=ji-1;
	166	jjp1=jj+1; jjm1=jj-1;
	167	summsk=(zsmask0(jip1,jj)+zsmask0(jim1,jj)+zsmask0(ji,jjp1)+zsmask0(ji,jjm1))
[5945]	168	IF (zdsmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
[5820]	169	sshn(ji,jj)=( zssh0(jip1,jj)*zsmask0(jip1,jj) &
	170	& + zssh0(jim1,jj)*zsmask0(jim1,jj) &
	171	& + zssh0(ji,jjp1)*zsmask0(ji,jjp1) &
	172	& + zssh0(ji,jjm1)*zsmask0(ji,jjm1))/summsk
[5823]	173	zsmask1(ji,jj)=1._wp
[9019]	174	ENDIF
[5790]	175	END DO
	176	END DO
[9019]	177	CALL lbc_lnk( sshn , 'T', 1._wp )
	178	CALL lbc_lnk( zsmask1, 'T', 1._wp )
[5790]	179	zssh0 = sshn
	180	zsmask0 = zsmask1
	181	END DO
	182	sshn(:,:) = sshn(:,:) * ssmask(:,:)
	183
	184	!=============================================================================
[6075]	185	!PM: Is this needed since introduction of VVL by default?
[9019]	186	IF ( .NOT.ln_linssh ) THEN
[5820]	187	! Reconstruction of all vertical scale factors at now time steps
	188	! =============================================================================
	189	! Horizontal scale factor interpolations
	190	! --------------------------------------
[5790]	191	DO jk = 1,jpk
[5820]	192	DO jj=1,jpj
	193	DO ji=1,jpi
	194	IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN
[8329]	195	e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) / &
	196	& ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) )
[5820]	197	ENDIF
	198	END DO
	199	END DO
[5790]	200	END DO
[9019]	201	!
[6074]	202	CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
	203	CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
	204	CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
[5790]	205
[9019]	206	! Vertical scale factor interpolations
	207	! ------------------------------------
[6074]	208	CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' )
	209	CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
	210	CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
[9019]	211
	212	! t- and w- points depth
	213	! ----------------------
[6074]	214	gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
	215	gdepw_n(:,:,1) = 0.0_wp
	216	gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
[5790]	217	DO jj = 1,jpj
	218	DO ji = 1,jpi
	219	DO jk = 2,mikt(ji,jj)-1
[6074]	220	gdept_n(ji,jj,jk) = gdept_0(ji,jj,jk)
	221	gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
	222	gde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj)
[5790]	223	END DO
[5945]	224	IF (mikt(ji,jj) > 1) THEN
[5790]	225	jk = mikt(ji,jj)
[6074]	226	gdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * e3w_n(ji,jj,jk)
	227	gdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk)
	228	gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk ) - sshn (ji,jj)
[5790]	229	END IF
	230	DO jk = mikt(ji,jj)+1, jpk
[6074]	231	gdept_n(ji,jj,jk) = gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)
	232	gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
	233	gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk ) - sshn (ji,jj)
[5790]	234	END DO
	235	END DO
	236	END DO
	237
[5823]	238	! t-, u- and v- water column thickness
	239	! ------------------------------------
[6075]	240	ht_n(:,:) = 0._wp ; hu_n(:,:) = 0._wp ; hv_n(:,:) = 0._wp
[5790]	241	DO jk = 1, jpkm1
[6074]	242	hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
	243	hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
	244	ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
[5790]	245	END DO
	246	! ! Inverse of the local depth
[6074]	247	r1_hu_n(:,:) = 1._wp / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)
	248	r1_hv_n(:,:) = 1._wp / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)
[5790]	249
	250	END IF
	251
	252	!=============================================================================
	253	! compute velocity
	254	! compute velocity in order to conserve barotropic velocity (modification by poderation of the scale factor).
	255	ub(:,:,:)=un(:,:,:)
	256	vb(:,:,:)=vn(:,:,:)
	257	DO jk = 1,jpk
[5920]	258	DO jj = 1,jpj
	259	DO ji = 1,jpi
[6074]	260	un(ji,jj,jk) = ub(ji,jj,jk)pe3u_b(ji,jj,jk)pumask_b(ji,jj,jk)/e3u_n(ji,jj,jk)*umask(ji,jj,jk);
	261	vn(ji,jj,jk) = vb(ji,jj,jk)pe3v_b(ji,jj,jk)pvmask_b(ji,jj,jk)/e3v_n(ji,jj,jk)*vmask(ji,jj,jk);
[5790]	262	END DO
	263	END DO
	264	END DO
[5802]	265
[5790]	266	! compute new velocity if we close a cell (check barotropic velocity and change velocity over the water column)
	267	! compute barotropic velocity now and after
	268	ztrp(:,:,:) = ub(:,:,:)*pe3u_b(:,:,:);
	269	zbub(:,:) = SUM(ztrp,DIM=3)
	270	ztrp(:,:,:) = vb(:,:,:)*pe3v_b(:,:,:);
	271	zbvb(:,:) = SUM(ztrp,DIM=3)
[6074]	272	ztrp(:,:,:) = un(:,:,:)*e3u_n(:,:,:);
[5790]	273	zbun(:,:) = SUM(ztrp,DIM=3)
[6074]	274	ztrp(:,:,:) = vn(:,:,:)*e3v_n(:,:,:);
[5790]	275	zbvn(:,:) = SUM(ztrp,DIM=3)
[5802]	276
[5820]	277	! new water column
[5945]	278	zhu1=0.0_wp ;
	279	zhv1=0.0_wp ;
[5790]	280	DO jk = 1,jpk
[6074]	281	zhu1(:,:) = zhu1(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
	282	zhv1(:,:) = zhv1(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
[5790]	283	END DO
[5820]	284
	285	! compute correction
[5790]	286	zucorr = 0._wp
	287	zvcorr = 0._wp
[5920]	288	DO jj = 1,jpj
	289	DO ji = 1,jpi
[5945]	290	IF (zbun(ji,jj) /= zbub(ji,jj) .AND. zhu1(ji,jj) /= 0._wp ) THEN
	291	zucorr(ji,jj) = (zbun(ji,jj) - zbub(ji,jj))/zhu1(ji,jj)
[5790]	292	END IF
[5945]	293	IF (zbvn(ji,jj) /= zbvb(ji,jj) .AND. zhv1(ji,jj) /= 0._wp ) THEN
	294	zvcorr(ji,jj) = (zbvn(ji,jj) - zbvb(ji,jj))/zhv1(ji,jj)
[5790]	295	END IF
	296	END DO
	297	END DO
[5820]	298
	299	! update velocity
[5790]	300	DO jk = 1,jpk
	301	un(:,:,jk)=(un(:,:,jk) - zucorr(:,:))*umask(:,:,jk)
	302	vn(:,:,jk)=(vn(:,:,jk) - zvcorr(:,:))*vmask(:,:,jk)
	303	END DO
	304
	305	!=============================================================================
	306	! compute temp and salt
	307	! compute new tn and sn if we open a new cell
	308	tsb (:,:,:,:) = tsn(:,:,:,:)
	309	zts0(:,:,:,:) = tsn(:,:,:,:)
	310	ztmask1(:,:,:) = ptmask_b(:,:,:)
	311	ztmask0(:,:,:) = ptmask_b(:,:,:)
[5920]	312	DO iz = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
[5790]	313	DO jk = 1,jpk-1
	314	zdmask=tmask(:,:,jk)-ztmask0(:,:,jk);
[5920]	315	DO jj = 2,jpj-1
	316	DO ji = fs_2,fs_jpim1
[5790]	317	jip1=ji+1; jim1=ji-1;
	318	jjp1=jj+1; jjm1=jj-1;
[5820]	319	summsk= (ztmask0(jip1,jj ,jk)+ztmask0(jim1,jj ,jk)+ztmask0(ji ,jjp1,jk)+ztmask0(ji ,jjm1,jk))
[5945]	320	IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp) THEN
[5820]	321	!! horizontal basic extrapolation
	322	tsn(ji,jj,jk,1)=( zts0(jip1,jj ,jk,1)*ztmask0(jip1,jj ,jk) &
	323	& +zts0(jim1,jj ,jk,1)*ztmask0(jim1,jj ,jk) &
	324	& +zts0(ji ,jjp1,jk,1)*ztmask0(ji ,jjp1,jk) &
	325	& +zts0(ji ,jjm1,jk,1)*ztmask0(ji ,jjm1,jk) ) / summsk
	326	tsn(ji,jj,jk,2)=( zts0(jip1,jj ,jk,2)*ztmask0(jip1,jj ,jk) &
	327	& +zts0(jim1,jj ,jk,2)*ztmask0(jim1,jj ,jk) &
	328	& +zts0(ji ,jjp1,jk,2)*ztmask0(ji ,jjp1,jk) &
	329	& +zts0(ji ,jjm1,jk,2)*ztmask0(ji ,jjm1,jk) ) / summsk
[5790]	330	ztmask1(ji,jj,jk)=1
[5945]	331	ELSEIF (zdmask(ji,jj) == 1._wp .AND. summsk == 0._wp) THEN
[5920]	332	!! vertical extrapolation if horizontal extrapolation failed
[5790]	333	jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
	334	summsk=(ztmask0(ji,jj,jkm1)+ztmask0(ji,jj,jkp1))
[5945]	335	IF (zdmask(ji,jj) == 1._wp .AND. summsk /= 0._wp ) THEN
[5790]	336	tsn(ji,jj,jk,1)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1) &
	337	& +zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1))/summsk
	338	tsn(ji,jj,jk,2)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1) &
	339	& +zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1))/summsk
[5823]	340	ztmask1(ji,jj,jk)=1._wp
[5790]	341	END IF
	342	END IF
	343	END DO
	344	END DO
	345	END DO
[5823]	346
	347	CALL lbc_lnk(tsn(:,:,:,1),'T',1._wp)
	348	CALL lbc_lnk(tsn(:,:,:,2),'T',1._wp)
[5945]	349	CALL lbc_lnk(ztmask1, 'T',1._wp)
[5823]	350
	351	! update
[5790]	352	zts0(:,:,:,:) = tsn(:,:,:,:)
	353	ztmask0 = ztmask1
[5823]	354
[5790]	355	END DO
[5823]	356
	357	! mask new tsn field
[5790]	358	tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)
	359	tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)
	360
[5802]	361	! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
[6075]	362	!PM: Is this IF needed since change to VVL by default
	363	IF (.NOT.ln_linssh) THEN
[5790]	364	DO jk = 2,jpk-1
	365	DO jj = 1,jpj
	366	DO ji = 1,jpi
[8329]	367	IF (zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND. &
	368	& (tmask(ji,jj,1)==0._wp .OR. ptmask_b(ji,jj,1)==0._wp) ) THEN
[5790]	369	!compute weight
[6074]	370	zdzp1 = MAX(0._wp,gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
	371	zdz = gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk )
	372	zdzm1 = MAX(0._wp,pdepw_b(ji,jj,jk ) - gdepw_n(ji,jj,jk ))
[5823]	373	IF (zdz .LT. 0._wp) THEN
[5820]	374	CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
[5790]	375	END IF
	376	tsn(ji,jj,jk,jp_tem) = ( zdzp1*tsb(ji,jj,jk+1,jp_tem) &
	377	& + zdz *tsb(ji,jj,jk ,jp_tem) &
[6074]	378	& + zdzm1*tsb(ji,jj,jk-1,jp_tem) )/e3t_n(ji,jj,jk)
[5790]	379	tsn(ji,jj,jk,jp_sal) = ( zdzp1*tsb(ji,jj,jk+1,jp_sal) &
	380	& + zdz *tsb(ji,jj,jk ,jp_sal) &
[6074]	381	& + zdzm1*tsb(ji,jj,jk-1,jp_sal) )/e3t_n(ji,jj,jk)
[5790]	382	END IF
	383	END DO
	384	END DO
	385	END DO
	386	END IF
	387
[5823]	388	! closed pool
	389	! -----------------------------------------------------------------------------------------
	390	! case we open a cell but no neigbour cells available to get an estimate of T and S
	391	WHERE (tmask(:,:,:) == 1._wp .AND. tsn(:,:,:,2) == 0._wp)
[5945]	392	tsn(:,:,:,2) = -99._wp ! Special value for closed pool (checking purpose in output.init)
[5823]	393	tmask(:,:,:) = 0._wp ! set mask to 0 to run
	394	umask(:,:,:) = 0._wp
	395	vmask(:,:,:) = 0._wp
[5790]	396	END WHERE
[5823]	397
	398	! set mbkt and mikt to 1 in thiese location
[5790]	399	WHERE (SUM(tmask,dim=3) == 0)
	400	mbkt(:,:)=1 ; mbku(:,:)=1 ; mbkv(:,:)=1
	401	mikt(:,:)=1 ; miku(:,:)=1 ; mikv(:,:)=1
	402	END WHERE
[5823]	403	! -------------------------------------------------------------------------------------------
[5790]	404	! compute new tn and sn if we close cell
	405	! nothing to do
[5802]	406	!
[5790]	407	END SUBROUTINE iscpl_rst_interpol
	408
[7646]	409	!!======================================================================
[5790]	410	END MODULE iscplrst

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

Download in other formats: