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/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90 @ 6069

Last change on this file since 6069 was 6069, checked in by timgraham, 9 years ago
Merge of dev_MetOffice_merge_2015 into branch (only NEMO directory for now).
File size: 20.7 KB

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

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

Download in other formats: