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.

isfcpl.F90 in NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/ISF – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/ISF/isfcpl.F90 @ 14644

Last change on this file since 14644 was 14644, checked in by sparonuz, 3 years ago
Merge trunk -r14642:HEAD
File size: 35.9 KB

Rev	Line
[11553]	1	MODULE isfcpl
	2	!!======================================================================
	3	!! * MODULE isfcpl *
	4	!!
	5	!! iceshelf coupling module : module managing the coupling between NEMO and an ice sheet model
	6	!!
	7	!!======================================================================
	8	!! History : 4.1 ! 2019-07 (P. Mathiot) Original code
	9	!!----------------------------------------------------------------------
	10
	11	!!----------------------------------------------------------------------
[14053]	12	!! isfrst : read/write iceshelf variables in/from restart
[11553]	13	!!----------------------------------------------------------------------
[14053]	14	USE oce ! ocean dynamics and tracers
	15	#if defined key_qco
	16	USE domqco , ONLY : dom_qco_zgr ! vertical scale factor interpolation
[14200]	17	#elif defined key_linssh
	18	! ! fix in time coordinate
[13237]	19	#else
[14053]	20	USE domvvl , ONLY : dom_vvl_zgr ! vertical scale factor interpolation
[13237]	21	#endif
[14053]	22	USE domutl , ONLY : dom_ngb ! find the closest grid point from a given lon/lat position
	23	USE isf_oce ! ice shelf variable
	24	USE isfutils, ONLY : debug
[11852]	25	!
[11553]	26	USE in_out_manager ! I/O manager
	27	USE iom ! I/O library
[14053]	28	USE lib_mpp , ONLY : mpp_sum, mpp_max ! mpp routine
[11553]	29	!
	30	IMPLICIT NONE
	31
	32	PRIVATE
	33
[14072]	34	PUBLIC isfcpl_rst_write, isfcpl_init ! iceshelf restart read and write
	35	PUBLIC isfcpl_ssh, isfcpl_tra, isfcpl_vol, isfcpl_cons ! iceshelf correction for ssh, tra, dyn and conservation
[11553]	36
	37	TYPE isfcons
[14053]	38	INTEGER :: ii ! i global
	39	INTEGER :: jj ! j global
	40	INTEGER :: kk ! k level
	41	REAL(wp):: dvol ! volume increment
	42	REAL(wp):: dsal ! salt increment
	43	REAL(wp):: dtem ! heat increment
	44	REAL(wp):: lon ! lon
	45	REAL(wp):: lat ! lat
	46	INTEGER :: ngb ! 0/1 (valid location or not (ie on halo or no neigbourg))
[11553]	47	END TYPE
	48	!
[12340]	49	!! * Substitutions
	50	# include "do_loop_substitute.h90"
[13237]	51	# include "domzgr_substitute.h90"
[14219]	52	# include "single_precision_substitute.h90"
	53
[11553]	54	!!----------------------------------------------------------------------
	55	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
	56	!! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
	57	!! Software governed by the CeCILL license (see ./LICENSE)
	58	!!----------------------------------------------------------------------
	59	CONTAINS
[12068]	60	SUBROUTINE isfcpl_init(Kbb, Kmm, Kaa)
[11902]	61	!!---------------------------------------------------------------------
[12077]	62	!! * ROUTINE iscpl_init *
[14072]	63	!!
	64	!! ** Purpose : correct ocean state for new wet cell and horizontal divergence
[12077]	65	!! correction for the dynamical adjustement
[11902]	66	!!
[12077]	67	!! ** Action : - compute ssh on new wet cell
	68	!! - compute T/S on new wet cell
	69	!! - compute horizontal divergence correction as a volume flux
	70	!! - compute the T/S/vol correction increment to keep trend to 0
	71	!!
[11902]	72	!!---------------------------------------------------------------------
[12068]	73	INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
	74	!!---------------------------------------------------------------------
[11902]	75	INTEGER :: id
	76	!!----------------------------------------------------------------------
[11852]	77	!
	78	! start on an euler time step
[12489]	79	l_1st_euler = .TRUE.
[14072]	80	!
[11902]	81	! allocation and initialisation to 0
[11553]	82	CALL isf_alloc_cpl()
	83	!
[11902]	84	! check presence of variable needed for coupling
	85	! iom_varid return 0 if not found
[11908]	86	id = 1
[11902]	87	id = id * iom_varid(numror, 'ssmask', ldstop = .false.)
	88	id = id * iom_varid(numror, 'tmask' , ldstop = .false.)
	89	id = id * iom_varid(numror, 'e3t_n' , ldstop = .false.)
	90	id = id * iom_varid(numror, 'e3u_n' , ldstop = .false.)
	91	id = id * iom_varid(numror, 'e3v_n' , ldstop = .false.)
	92	IF(lwp) WRITE(numout,*) ' isfcpl_init:', id
	93	IF (id == 0) THEN
[14219]	94	IF(lwp) WRITE(numout,*) ' isfcpl_init: restart variables for ice sheet coupling are missing, skip coupling for this leg '
[11902]	95	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
	96	IF(lwp) WRITE(numout,*) ''
	97	ELSE
	98	! extrapolation ssh
[12068]	99	CALL isfcpl_ssh(Kbb, Kmm, Kaa)
[11902]	100	!
	101	! extrapolation tracer properties
[12068]	102	CALL isfcpl_tra(Kmm)
[11902]	103	!
	104	! correction of the horizontal divergence and associated temp. and salt content flux
	105	! Need to : - include in the cpl cons the risfcpl_vol/tsc contribution
	106	! - decide how to manage thickness level change in conservation
[12068]	107	CALL isfcpl_vol(Kmm)
[11902]	108	!
	109	! apply the 'conservation' method
[12068]	110	IF ( ln_isfcpl_cons ) CALL isfcpl_cons(Kmm)
[11902]	111	!
	112	END IF
[11553]	113	!
	114	! mask velocity properly (mask used in restart not compatible with new mask)
[12068]	115	uu(:,:,:,Kmm) = uu(:,:,:,Kmm) * umask(:,:,:)
	116	vv(:,:,:,Kmm) = vv(:,:,:,Kmm) * vmask(:,:,:)
[11553]	117	!
[11852]	118	! all before fields set to now values
[12068]	119	ts (:,:,:,:,Kbb) = ts (:,:,:,:,Kmm)
	120	uu (:,:,:,Kbb) = uu (:,:,:,Kmm)
	121	vv (:,:,:,Kbb) = vv (:,:,:,Kmm)
	122	ssh (:,:,Kbb) = ssh (:,:,Kmm)
[14200]	123	#if ! defined key_qco && ! defined key_linssh
[12068]	124	e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
[14072]	125	#endif
[11553]	126	END SUBROUTINE isfcpl_init
[14053]	127
[14072]	128
[14053]	129	SUBROUTINE isfcpl_rst_write( kt, Kmm )
[11553]	130	!!---------------------------------------------------------------------
[11931]	131	!! * ROUTINE iscpl_rst_write *
[14072]	132	!!
[11931]	133	!! ** Purpose : write icesheet coupling variables in restart
	134	!!
[11553]	135	!!-------------------------- IN --------------------------------------
	136	INTEGER, INTENT(in) :: kt
[12068]	137	INTEGER, INTENT(in) :: Kmm ! ocean time level index
[11553]	138	!!----------------------------------------------------------------------
[13237]	139	INTEGER :: jk ! loop index
[14053]	140	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, ze3u, ze3v, zgdepw ! for qco substitution
[13237]	141	!!----------------------------------------------------------------------
[11553]	142	!
[13237]	143	DO jk = 1, jpk
	144	ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
	145	ze3u(:,:,jk) = e3u(:,:,jk,Kmm)
	146	ze3v(:,:,jk) = e3v(:,:,jk,Kmm)
	147	!
	148	zgdepw(:,:,jk) = gdepw(:,:,jk,Kmm)
[14072]	149	END DO
[13237]	150	!
[13970]	151	CALL iom_rstput( kt, nitrst, numrow, 'tmask' , tmask )
	152	CALL iom_rstput( kt, nitrst, numrow, 'ssmask' , ssmask )
	153	CALL iom_rstput( kt, nitrst, numrow, 'e3t_n' , ze3t )
	154	CALL iom_rstput( kt, nitrst, numrow, 'e3u_n' , ze3u )
	155	CALL iom_rstput( kt, nitrst, numrow, 'e3v_n' , ze3v )
	156	CALL iom_rstput( kt, nitrst, numrow, 'gdepw_n', zgdepw )
[11553]	157	!
	158	END SUBROUTINE isfcpl_rst_write
	159
[14072]	160
[12068]	161	SUBROUTINE isfcpl_ssh(Kbb, Kmm, Kaa)
[14072]	162	!!----------------------------------------------------------------------
[11553]	163	!! * ROUTINE iscpl_ssh *
[14072]	164	!!
[11931]	165	!! ** Purpose : basic guess of ssh in new wet cell
[14072]	166	!!
[11553]	167	!! ** Method : basic extrapolation from neigbourg cells
	168	!!
	169	!!----------------------------------------------------------------------
	170	!!
[12068]	171	INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
	172	!!----------------------------------------------------------------------
[11908]	173	INTEGER :: ji, jj, jd, jk !! loop index
[11553]	174	INTEGER :: jip1, jim1, jjp1, jjm1
	175	!!
	176	REAL(wp):: zsummsk
	177	REAL(wp), DIMENSION(jpi,jpj) :: zdssmask, zssmask0, zssmask_b, zssh
	178	!!----------------------------------------------------------------------
	179	!
[13970]	180	CALL iom_get( numror, jpdom_auto, 'ssmask' , zssmask_b ) ! need to extrapolate T/S
[11553]	181
[14072]	182	! compute new ssh if we open a full water column
[11553]	183	! rude average of the closest neigbourgs (e1e2t not taking into account)
	184	!
[12068]	185	zssh(:,:) = ssh(:,:,Kmm)
[11553]	186	zssmask0(:,:) = zssmask_b(:,:)
	187	!
	188	DO jd = 1, nn_drown
	189	!
	190	zdssmask(:,:) = ssmask(:,:) - zssmask0(:,:)
[13295]	191	DO_2D( 0, 0, 0, 0 )
[14053]	192	jip1=ji+1 ; jim1=ji-1
	193	jjp1=jj+1 ; jjm1=jj-1
[12353]	194	!
	195	zsummsk = zssmask0(jip1,jj) + zssmask0(jim1,jj) + zssmask0(ji,jjp1) + zssmask0(ji,jjm1)
	196	!
	197	IF (zdssmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp) THEN
	198	ssh(ji,jj,Kmm)=( zssh(jip1,jj)*zssmask0(jip1,jj) &
[14053]	199	& + zssh(jim1,jj)*zssmask0(jim1,jj) &
	200	& + zssh(ji,jjp1)*zssmask0(ji,jjp1) &
	201	& + zssh(ji,jjm1)*zssmask0(ji,jjm1)) / zsummsk
[12353]	202	zssmask_b(ji,jj) = 1._wp
	203	ENDIF
	204	END_2D
[11553]	205	!
[12068]	206	zssh(:,:) = ssh(:,:,Kmm)
[11553]	207	zssmask0(:,:) = zssmask_b(:,:)
	208	!
[14644]	209	CALL lbc_lnk( 'iscplrst', zssh, 'T', 1.0_wp, zssmask0, 'T', 1.0_wp )
[11553]	210	!
	211	END DO
	212	!
[12068]	213	! update ssh(:,:,Kmm)
	214	ssh(:,:,Kmm) = zssh(:,:) * ssmask(:,:)
[11553]	215	!
[12068]	216	ssh(:,:,Kbb) = ssh(:,:,Kmm)
[11553]	217	!
[14219]	218	IF ( ln_isfdebug ) CALL debug('isfcpl_ssh: sshn',CASTWP(ssh(:,:,Kmm)))
[11902]	219	!
[11823]	220	! recompute the vertical scale factor, depth and water thickness
[12068]	221	IF(lwp) write(numout,*) 'isfcpl_ssh : recompute scale factor from ssh (new wet cell,Kmm)'
[11908]	222	IF(lwp) write(numout,*) '~~~~~~~~~~~'
[14200]	223	#if defined key_qco
	224	CALL dom_qco_zgr(Kbb, Kmm)
	225	#elif defined key_linssh
	226	! linear ssh : fix in time coord.
	227	#else
[11908]	228	DO jk = 1, jpk
[13237]	229	e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( 1._wp + (ht_0(:,:) + ssh(:,:,Kmm)) * r1_ht_0(:,:) )
[11908]	230	END DO
[12068]	231	e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
	232	CALL dom_vvl_zgr(Kbb, Kmm, Kaa)
[13237]	233	#endif
[11823]	234	!
[11553]	235	END SUBROUTINE isfcpl_ssh
	236
[14072]	237
[12068]	238	SUBROUTINE isfcpl_tra(Kmm)
[14072]	239	!!----------------------------------------------------------------------
[11553]	240	!! * ROUTINE iscpl_tra *
[14072]	241	!!
	242	!! ** Purpose : compute new tn, sn in case of evolving geometry of ice shelves
	243	!!
[11553]	244	!! ** Method : tn, sn : basic extrapolation from neigbourg cells
	245	!!
	246	!!----------------------------------------------------------------------
[12068]	247	INTEGER, INTENT(in) :: Kmm ! ocean time level index
	248	!!----------------------------------------------------------------------
[11553]	249	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b
	250	!REAL(wp), DIMENSION(:,:,: ), INTENT(in ) :: pdepw_b !! depth w before
	251	!!
	252	INTEGER :: ji, jj, jk, jd !! loop index
	253	INTEGER :: jip1, jim1, jjp1, jjm1, jkp1, jkm1
	254	!!
	255	REAL(wp):: zsummsk
	256	REAL(wp):: zdz, zdzm1, zdzp1
	257	!!
[14072]	258	REAL(wp), DIMENSION(jpi,jpj) :: zdmask
[11553]	259	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask0, zwmaskn
	260	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask1, zwmaskb, ztmp3d
	261	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) :: zts0
	262	!!----------------------------------------------------------------------
[14072]	263	!
[13970]	264	CALL iom_get( numror, jpdom_auto, 'tmask' , ztmask_b ) ! need to extrapolate T/S
	265	!CALL iom_get( numror, jpdom_auto, 'wmask' , zwmask_b ) ! need to extrapolate T/S
	266	!CALL iom_get( numror, jpdom_auto, 'gdepw_n', zdepw_b(:,:,:) ) ! need to interpol vertical profile (vvl)
[11553]	267	!
[14072]	268	!
[11553]	269	! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
	270	!PM: Is this IF needed since change to VVL by default
	271	!bugged : to be corrected (PM)
	272	! back up original t/s/mask
[12068]	273	!tsb (:,:,:,:) = ts(:,:,:,:,Kmm)
[11553]	274	!
	275	! compute new T/S (interpolation) if vvl only for common wet cell in before and after wmask
	276
	277	! IF (.NOT.ln_linssh) THEN
	278	! DO jk = 2,jpk-1
	279	! DO jj = 1,jpj
	280	! DO ji = 1,jpi
	281	! IF (wmask(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND. (tmask(ji,jj,1)==0._wp .OR. ztmask_b(ji,jj,1)==0._wp) ) THEN
	282	!
	283	! !compute weight
[12068]	284	! zdzp1 = MAX(0._wp,pdepw_b(ji,jj,jk+1) - gdepw(ji,jj,jk+1,Kmm))
	285	! zdzm1 = MAX(0._wp,gdepw(ji,jj,jk ,Kmm) - pdepw_b(ji,jj,jk ))
	286	! zdz = e3t(ji,jj,jk,Kmm) - zdzp1 - zdzm1 ! if isf : e3t = gdepw(ji,jj,jk+1,Kmm)- gdepw(ji,jj,jk,Kmm)
[11553]	287	!
	288	! IF (zdz .LT. 0._wp) THEN
	289	! CALL ctl_stop( 'STOP', 'rst_iscpl : unable to compute the interpolation' )
	290	! END IF
	291	!
[12068]	292	! ts(ji,jj,jk,jp_tem,Kmm) = ( zdzp1*ts(ji,jj,jk+1,jp_tem,Kbb) &
	293	! & + zdz *ts(ji,jj,jk ,jp_tem,Kbb) &
	294	! & + zdzm1*ts(ji,jj,jk-1,jp_tem,Kbb) )/e3t(ji,jj,jk,Kmm)
[11553]	295	!
[12068]	296	! ts(ji,jj,jk,jp_sal,Kmm) = ( zdzp1*ts(ji,jj,jk+1,jp_sal,Kbb) &
	297	! & + zdz *ts(ji,jj,jk ,jp_sal,Kbb) &
	298	! & + zdzm1*ts(ji,jj,jk-1,jp_sal,Kbb) )/e3t(ji,jj,jk,Kmm)
[11553]	299	!
	300	! END IF
	301	! END DO
	302	! END DO
	303	! END DO
	304	! END IF
	305
[12068]	306	zts0(:,:,:,:) = ts(:,:,:,:,Kmm)
[11553]	307	ztmask0(:,:,:) = ztmask_b(:,:,:)
	308	ztmask1(:,:,:) = ztmask_b(:,:,:)
	309	!
	310	! iterate the extrapolation processes nn_drown times
	311	DO jd = 1,nn_drown ! resolution dependent (OK for ISOMIP+ case)
	312	DO jk = 1,jpk-1
	313	!
	314	! define new wet cell
	315	zdmask(:,:) = tmask(:,:,jk) - ztmask0(:,:,jk);
	316	!
[13295]	317	DO_2D( 0, 0, 0, 0 )
[12353]	318	jip1=ji+1; jim1=ji-1;
	319	jjp1=jj+1; jjm1=jj-1;
	320	!
	321	! check if a wet neigbourg cell is present
	322	zsummsk = ztmask0(jip1,jj ,jk) + ztmask0(jim1,jj ,jk) &
	323	+ ztmask0(ji ,jjp1,jk) + ztmask0(ji ,jjm1,jk)
	324	!
	325	! if neigbourg wet cell available at the same level
	326	IF ( zdmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp ) THEN
[11553]	327	!
[12353]	328	! horizontal basic extrapolation
	329	ts(ji,jj,jk,1,Kmm)=( zts0(jip1,jj ,jk,1) * ztmask0(jip1,jj ,jk) &
	330	& + zts0(jim1,jj ,jk,1) * ztmask0(jim1,jj ,jk) &
	331	& + zts0(ji ,jjp1,jk,1) * ztmask0(ji ,jjp1,jk) &
	332	& + zts0(ji ,jjm1,jk,1) * ztmask0(ji ,jjm1,jk) ) / zsummsk
	333	ts(ji,jj,jk,2,Kmm)=( zts0(jip1,jj ,jk,2) * ztmask0(jip1,jj ,jk) &
	334	& + zts0(jim1,jj ,jk,2) * ztmask0(jim1,jj ,jk) &
	335	& + zts0(ji ,jjp1,jk,2) * ztmask0(ji ,jjp1,jk) &
	336	& + zts0(ji ,jjm1,jk,2) * ztmask0(ji ,jjm1,jk) ) / zsummsk
	337	!
	338	! update mask for next pass
	339	ztmask1(ji,jj,jk)=1
	340	!
	341	! in case no neigbourg wet cell available at the same level
	342	! check if a wet cell is available below
	343	ELSEIF (zdmask(ji,jj) == 1._wp .AND. zsummsk == 0._wp) THEN
	344	!
	345	! vertical extrapolation if horizontal extrapolation failed
	346	jkm1=max(1,jk-1) ; jkp1=min(jpk,jk+1)
	347	!
[11553]	348	! check if a wet neigbourg cell is present
[12353]	349	zsummsk = ztmask0(ji,jj,jkm1) + ztmask0(ji,jj,jkp1)
	350	IF (zdmask(ji,jj) == 1._wp .AND. zsummsk /= 0._wp ) THEN
	351	ts(ji,jj,jk,1,Kmm)=( zts0(ji,jj,jkp1,1)*ztmask0(ji,jj,jkp1) &
	352	& + zts0(ji,jj,jkm1,1)*ztmask0(ji,jj,jkm1)) / zsummsk
	353	ts(ji,jj,jk,2,Kmm)=( zts0(ji,jj,jkp1,2)*ztmask0(ji,jj,jkp1) &
	354	& + zts0(ji,jj,jkm1,2)*ztmask0(ji,jj,jkm1)) / zsummsk
[11553]	355	!
	356	! update mask for next pass
[12353]	357	ztmask1(ji,jj,jk)=1._wp
[11553]	358	END IF
[12353]	359	END IF
	360	END_2D
[11553]	361	END DO
	362	!
	363	! update temperature and salinity and mask
[12068]	364	zts0(:,:,:,:) = ts(:,:,:,:,Kmm)
[11553]	365	ztmask0(:,:,:) = ztmask1(:,:,:)
	366	!
[14644]	367	CALL lbc_lnk( 'iscplrst', zts0(:,:,:,jp_tem), 'T', 1.0_wp, zts0(:,:,:,jp_sal), 'T', 1.0_wp, ztmask0, 'T', 1.0_wp)
[11553]	368	!
	369	END DO ! nn_drown
	370	!
[12068]	371	! mask new ts(:,:,:,:,Kmm) field
	372	ts(:,:,:,jp_tem,Kmm) = zts0(:,:,:,jp_tem) * tmask(:,:,:)
	373	ts(:,:,:,jp_sal,Kmm) = zts0(:,:,:,jp_sal) * tmask(:,:,:)
[11553]	374	!
	375	! sanity check
	376	! -----------------------------------------------------------------------------------------
	377	! case we open a cell but no neigbour cells available to get an estimate of T and S
[13295]	378	DO_3D( 1, 1, 1, 1, 1,jpk-1 )
[12340]	379	IF (tmask(ji,jj,jk) == 1._wp .AND. ts(ji,jj,jk,2,Kmm) == 0._wp) &
	380	& CALL ctl_stop('STOP', 'failing to fill all new weet cell, &
	381	& try increase nn_drown or activate XXXX &
	382	& in your domain cfg computation' )
	383	END_3D
[14072]	384	!
[11553]	385	END SUBROUTINE isfcpl_tra
	386
[14072]	387
[12068]	388	SUBROUTINE isfcpl_vol(Kmm)
[14072]	389	!!----------------------------------------------------------------------
[11553]	390	!! * ROUTINE iscpl_vol *
[14072]	391	!!
	392	!! ** Purpose : compute the correction of the local divergence to apply
[11931]	393	!! during the first time step after the coupling.
	394	!!
	395	!! ** Method : - compute horizontal vol div. before/after coupling
	396	!! - compute vertical input
	397	!! - compute correction
[14072]	398	!!
[11553]	399	!!----------------------------------------------------------------------
	400	!!
[12068]	401	INTEGER, INTENT(in) :: Kmm ! ocean time level index
	402	!!----------------------------------------------------------------------
[14072]	403	INTEGER :: ji, jj, jk
[11553]	404	INTEGER :: ikb, ikt
	405	!!
[11931]	406	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zqvolb, zqvoln ! vol flux div. before/after coupling
	407	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3u_b, ze3v_b ! vertical scale factor before/after coupling
	408	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b ! mask before coupling
[11553]	409	!!----------------------------------------------------------------------
	410	!
[13970]	411	CALL iom_get( numror, jpdom_auto, 'tmask' , ztmask_b )
	412	CALL iom_get( numror, jpdom_auto, 'e3u_n' , ze3u_b )
	413	CALL iom_get( numror, jpdom_auto, 'e3v_n' , ze3v_b )
[11553]	414	!
[11931]	415	! 1.0: compute horizontal volume flux divergence difference before-after coupling
	416	!
[11553]	417	DO jk = 1, jpk ! Horizontal slab
[11931]	418	! 1.1: get volume flux before coupling (>0 out)
[13295]	419	DO_2D( 0, 0, 0, 0 )
[13237]	420	zqvolb(ji,jj,jk) = &
	421	& ( e2u(ji ,jj ) * ze3u_b(ji ,jj ,jk) * uu(ji ,jj ,jk,Kmm) &
	422	& - e2u(ji-1,jj ) * ze3u_b(ji-1,jj ,jk) * uu(ji-1,jj ,jk,Kmm) &
	423	& + e1v(ji ,jj ) * ze3v_b(ji ,jj ,jk) * vv(ji ,jj ,jk,Kmm) &
	424	& - e1v(ji ,jj-1) * ze3v_b(ji ,jj-1,jk) * vv(ji ,jj-1,jk,Kmm) ) &
[12340]	425	& * ztmask_b(ji,jj,jk)
	426	END_2D
[11553]	427	!
[11931]	428	! 1.2: get volume flux after coupling (>0 out)
[14072]	429	! properly mask velocity
[11553]	430	! (velocity are still mask with old mask at this stage)
[12068]	431	uu(:,:,jk,Kmm) = uu(:,:,jk,Kmm) * umask(:,:,jk)
	432	vv(:,:,jk,Kmm) = vv(:,:,jk,Kmm) * vmask(:,:,jk)
[11931]	433	! compute volume flux divergence after coupling
[13295]	434	DO_2D( 0, 0, 0, 0 )
[13237]	435	zqvoln(ji,jj,jk) = &
	436	& ( e2u(ji ,jj ) * e3u(ji ,jj ,jk,Kmm) * uu(ji ,jj ,jk,Kmm) &
	437	& - e2u(ji-1,jj ) * e3u(ji-1,jj ,jk,Kmm) * uu(ji-1,jj ,jk,Kmm) &
	438	& + e1v(ji ,jj ) * e3v(ji ,jj ,jk,Kmm) * vv(ji ,jj ,jk,Kmm) &
	439	& - e1v(ji ,jj-1) * e3v(ji ,jj-1,jk,Kmm) * vv(ji ,jj-1,jk,Kmm) ) &
[12340]	440	& * tmask(ji,jj,jk)
	441	END_2D
[11553]	442	!
[11931]	443	! 1.3: get 3d volume flux difference (before - after cpl) (>0 out)
	444	! correction to add is _b - _n
[11553]	445	risfcpl_vol(:,:,jk) = zqvolb(:,:,jk) - zqvoln(:,:,jk)
	446	END DO
	447	!
[11931]	448	! 2.0: include the contribution of the vertical velocity in the volume flux correction
	449	!
[13295]	450	DO_2D( 0, 0, 0, 0 )
[12340]	451	!
	452	ikt = mikt(ji,jj)
	453	IF ( ikt > 1 .AND. ssmask(ji,jj) == 1 ) THEN
	454	risfcpl_vol(ji,jj,ikt) = risfcpl_vol(ji,jj,ikt) + SUM(zqvolb(ji,jj,1:ikt-1)) ! test sign
	455	ENDIF
	456	!
	457	END_2D
[11553]	458	!
[13226]	459	CALL lbc_lnk( 'iscpl', risfcpl_vol, 'T', 1.0_wp )
[11931]	460	!
[12068]	461	! 3.0: set total correction (div, tr(:,:,:,:,Krhs), ssh)
[11931]	462	!
	463	! 3.1: mask volume flux divergence correction
[11553]	464	risfcpl_vol(:,:,:) = risfcpl_vol(:,:,:) * tmask(:,:,:)
	465	!
[12068]	466	! 3.2: get 3d tr(:,:,:,:,Krhs) increment to apply at the first time step
[14072]	467	! temperature and salt content flux computed using local ts(:,:,:,:,Kmm)
[11553]	468	! (very simple advection scheme)
	469	! (>0 out)
[12068]	470	risfcpl_tsc(:,:,:,jp_tem) = -risfcpl_vol(:,:,:) * ts(:,:,:,jp_tem,Kmm)
	471	risfcpl_tsc(:,:,:,jp_sal) = -risfcpl_vol(:,:,:) * ts(:,:,:,jp_sal,Kmm)
[11553]	472	!
[11931]	473	! 3.3: ssh correction (for dynspg_ts)
[11553]	474	risfcpl_ssh(:,:) = 0.0
	475	DO jk = 1,jpk
	476	risfcpl_ssh(:,:) = risfcpl_ssh(:,:) + risfcpl_vol(:,:,jk) * r1_e1e2t(:,:)
	477	END DO
[14053]	478	!
[11553]	479	END SUBROUTINE isfcpl_vol
	480
[14072]	481
[12068]	482	SUBROUTINE isfcpl_cons(Kmm)
[14072]	483	!!----------------------------------------------------------------------
[11553]	484	!! * ROUTINE iscpl_cons *
[14072]	485	!!
[11553]	486	!! ** Purpose : compute the corrective increment in volume/salt/heat to put back the vol/heat/salt
	487	!! removed or added during the coupling processes (wet or dry new cell)
[14072]	488	!!
[11553]	489	!! ** Method : - compare volume/heat/salt before and after
	490	!! - look for the closest wet cells (share amoung neigbourgs if there are)
	491	!! - build the correction increment to applied at each time step
[14072]	492	!!
[11553]	493	!!----------------------------------------------------------------------
	494	!
[11931]	495	TYPE(isfcons), DIMENSION(:),ALLOCATABLE :: zisfpts ! list of point receiving a correction
[11553]	496	!
[12068]	497	!!----------------------------------------------------------------------
	498	INTEGER, INTENT(in) :: Kmm ! ocean time level index
	499	!!----------------------------------------------------------------------
[11931]	500	INTEGER :: ji , jj , jk , jproc ! loop index
	501	INTEGER :: jip1 , jim1, jjp1, jjm1 ! dummy indices
	502	INTEGER :: iig , ijg, ik ! dummy indices
	503	INTEGER :: jisf ! start, end and current position in the increment array
[14072]	504	INTEGER :: ingb, ifind ! 0/1 target found or need to be found
	505	INTEGER :: nisfl_area ! global number of cell concerned by the wet->dry case
[11931]	506	INTEGER, DIMENSION(jpnij) :: nisfl ! local number of cell concerned by the wet->dry case
[11553]	507	!
	508	REAL(wp) :: z1_sum, z1_rdtiscpl
[11931]	509	REAL(wp) :: zdtem, zdsal, zdvol, zratio ! tem, sal, vol increment
[14072]	510	REAL(wp) :: zlon , zlat ! target location
[11931]	511	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztmask_b ! mask before
	512	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t_b ! scale factor before
	513	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zt_b ! scale factor before
	514	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zs_b ! scale factor before
[11553]	515	!!----------------------------------------------------------------------
	516
[11931]	517	!==============================================================================
	518	! 1.0: initialisation
	519	!==============================================================================
	520
[11553]	521	! get restart variable
[13970]	522	CALL iom_get( numror, jpdom_auto, 'tmask' , ztmask_b(:,:,:) ) ! need to extrapolate T/S
	523	CALL iom_get( numror, jpdom_auto, 'e3t_n' , ze3t_b(:,:,:) )
	524	CALL iom_get( numror, jpdom_auto, 'tn' , zt_b(:,:,:) )
	525	CALL iom_get( numror, jpdom_auto, 'sn' , zs_b(:,:,:) )
[11553]	526
	527	! compute run length
	528	nstp_iscpl = nitend - nit000 + 1
[12489]	529	rdt_iscpl = nstp_iscpl * rn_Dt
[14072]	530	z1_rdtiscpl = 1._wp / rdt_iscpl
[11931]	531
[11553]	532	IF (lwp) WRITE(numout,*) ' nb of stp for cons = ', nstp_iscpl
	533	IF (lwp) WRITE(numout,*) ' coupling time step = ', rdt_iscpl
	534
[11931]	535	! initialisation correction
	536	risfcpl_cons_vol = 0.0
	537	risfcpl_cons_ssh = 0.0
	538	risfcpl_cons_tsc = 0.0
[11553]	539
	540	!==============================================================================
[11931]	541	! 2.0: diagnose the heat, salt and volume input and compute the correction variable
	542	! for case where we wet a cell or cell still wet (no change in cell status)
[11553]	543	!==============================================================================
	544
	545	DO jk = 1,jpk-1
[13286]	546	DO jj = Njs0,Nje0
	547	DO ji = Nis0,Nie0
[11553]	548
	549	! volume diff
[13237]	550	zdvol = e3t (ji,jj,jk,Kmm) * tmask (ji,jj,jk) &
	551	& - ze3t_b(ji,jj,jk ) * ztmask_b(ji,jj,jk)
[11553]	552
	553	! heat diff
[12068]	554	zdtem = ts (ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm) * tmask (ji,jj,jk) &
[11931]	555	- zt_b(ji,jj,jk) * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
[11553]	556
	557	! salt diff
[12068]	558	zdsal = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm) * tmask (ji,jj,jk) &
[11931]	559	- zs_b(ji,jj,jk) * ze3t_b(ji,jj,jk) * ztmask_b(ji,jj,jk)
[14072]	560
[11553]	561	! volume, heat and salt differences in each cell (>0 means correction is an outward flux)
[11931]	562	! in addition to the geometry change unconservation, need to add the divergence correction as it is flux across the boundary
	563	risfcpl_cons_vol(ji,jj,jk) = ( zdvol * e1e2t(ji,jj) + risfcpl_vol(ji,jj,jk) ) * z1_rdtiscpl
	564	risfcpl_cons_tsc(ji,jj,jk,jp_sal) = ( - zdsal * e1e2t(ji,jj) + risfcpl_tsc(ji,jj,jk,jp_sal) ) * z1_rdtiscpl
	565	risfcpl_cons_tsc(ji,jj,jk,jp_tem) = ( - zdtem * e1e2t(ji,jj) + risfcpl_tsc(ji,jj,jk,jp_tem) ) * z1_rdtiscpl
[11553]	566
	567	END DO
	568	END DO
	569	END DO
	570	!
[11931]	571	!==============================================================================
	572	! 3.0: diagnose the heat, salt and volume input and compute the correction variable
	573	! for case where we close a cell
	574	!==============================================================================
	575	!
[11553]	576	! compute the total number of point receiving a correction increment for each processor
	577	! local
[11931]	578	nisfl(:)=0
[11553]	579	DO jk = 1,jpk-1
[13286]	580	DO jj = Njs0,Nje0
	581	DO ji = Nis0,Nie0
[11553]	582	jip1=MIN(ji+1,jpi) ; jim1=MAX(ji-1,1) ; jjp1=MIN(jj+1,jpj) ; jjm1=MAX(jj-1,1) ;
[14072]	583	IF ( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) THEN
[13237]	584	nisfl(narea) = nisfl(narea) + MAX(SUM(tmask(jim1:jip1,jjm1:jjp1,jk)),1._wp)
	585	ENDIF
[11553]	586	ENDDO
	587	ENDDO
	588	ENDDO
	589	!
[14286]	590	! global
[11931]	591	CALL mpp_sum('isfcpl',nisfl )
[11553]	592	!
	593	! allocate list of point receiving correction
[11931]	594	ALLOCATE(zisfpts(nisfl(narea)))
	595	!
[14219]	596	zisfpts(:) = isfcons(0,0,0,-HUGE(1.0_wp), -HUGE(1.0_wp), -HUGE(1.0_wp), -HUGE(1.0_wp), -HUGE(1.0_wp), 0)
[11553]	597	!
	598	! start computing the correction and fill zisfpts
	599	! local
[11931]	600	jisf = 0
[11553]	601	DO jk = 1,jpk-1
[13286]	602	DO jj = Njs0,Nje0
	603	DO ji = Nis0,Nie0
[11553]	604	IF ( tmask(ji,jj,jk) == 0._wp .AND. ztmask_b(ji,jj,jk) == 1._wp ) THEN
	605
	606	jip1=MIN(ji+1,jpi) ; jim1=MAX(ji-1,1) ; jjp1=MIN(jj+1,jpj) ; jjm1=MAX(jj-1,1) ;
	607
	608	zdvol = risfcpl_cons_vol(ji,jj,jk )
	609	zdsal = risfcpl_cons_tsc(ji,jj,jk,jp_sal)
	610	zdtem = risfcpl_cons_tsc(ji,jj,jk,jp_tem)
	611
	612	IF ( SUM( tmask(jim1:jip1,jjm1:jjp1,jk) ) > 0._wp ) THEN
	613	! spread correction amoung neigbourg wet cells (horizontal direction first)
	614	! as it is a rude correction corner and lateral cell have the same weight
	615	!
	616	z1_sum = 1._wp / SUM( tmask(jim1:jip1,jjm1:jjp1,jk) )
	617	!
	618	! lateral cells
	619	IF (tmask(jip1,jj ,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jj , jk, zdvol, zdsal, zdtem, z1_sum)
	620	IF (tmask(jim1,jj ,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jj , jk, zdvol, zdsal, zdtem, z1_sum)
	621	IF (tmask(ji ,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, ji , jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
	622	IF (tmask(ji ,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, ji , jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
	623	!
	624	! corner cells
	625	IF (tmask(jip1,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
	626	IF (tmask(jim1,jjm1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jjm1, jk, zdvol, zdsal, zdtem, z1_sum)
	627	IF (tmask(jim1,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jim1, jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
	628	IF (tmask(jip1,jjp1,jk) == 1) CALL update_isfpts(zisfpts, jisf, jip1, jjp1, jk, zdvol, zdsal, zdtem, z1_sum)
	629	!
	630	ELSE IF ( tmask(ji,jj,jk+1) == 1._wp ) THEN
	631	! spread correction amoung neigbourg wet cells (vertical direction)
[13226]	632	CALL update_isfpts(zisfpts, jisf, ji , jj , jk+1, zdvol, zdsal, zdtem, 1.0_wp, 0)
[11553]	633	ELSE
	634	! need to find where to put correction in later on
[13226]	635	CALL update_isfpts(zisfpts, jisf, ji , jj , jk , zdvol, zdsal, zdtem, 1.0_wp, 1)
[11553]	636	END IF
	637	END IF
	638	END DO
	639	END DO
	640	END DO
	641	!
	642	! share data among all processes because for some point we need to find the closest wet point (could be on other process)
[11931]	643	DO jproc=1,jpnij
[14072]	644	!
[11931]	645	! share total number of isf point treated for proc jproc
	646	IF (jproc==narea) THEN
	647	nisfl_area=nisfl(jproc)
	648	ELSE
	649	nisfl_area=0
	650	END IF
	651	CALL mpp_max('isfcpl',nisfl_area)
[11553]	652	!
[11931]	653	DO jisf = 1,nisfl_area
	654	!
	655	IF (jproc==narea) THEN
	656	! indices (conversion to global indices and sharing)
	657	iig = zisfpts(jisf)%ii ; ijg = zisfpts(jisf)%jj ; ik = zisfpts(jisf)%kk
	658	!
	659	! data
	660	zdvol = zisfpts(jisf)%dvol ; zdsal = zisfpts(jisf)%dsal ; zdtem = zisfpts(jisf)%dtem
	661	!
	662	! location
	663	zlat = zisfpts(jisf)%lat ; zlon = zisfpts(jisf)%lon
	664	!
	665	! find flag
	666	ingb = zisfpts(jisf)%ngb
	667	ELSE
[14072]	668	iig =0 ; ijg =0 ; ik =0
[14219]	669	zdvol=-HUGE(1.0_wp) ; zdsal=-HUGE(1.0_wp) ; zdtem=-HUGE(1.0)
	670	zlat =-HUGE(1.0_wp) ; zlon =-HUGE(1.0_wp)
[11931]	671	ingb = 0
	672	END IF
	673	!
	674	! share data (need synchronisation of data as get_correction call a global com)
	675	CALL mpp_max('isfcpl',iig) ; CALL mpp_max('isfcpl',ijg) ; CALL mpp_max('isfcpl',ik)
	676	CALL mpp_max('isfcpl',zdvol) ; CALL mpp_max('isfcpl',zdsal) ; CALL mpp_max('isfcpl',zdtem)
	677	CALL mpp_max('isfcpl',zlat) ; CALL mpp_max('isfcpl',zlon)
	678	CALL mpp_max('isfcpl',ingb)
	679	!
	680	! fill the 3d correction array
	681	CALL get_correction(iig, ijg, ik, zlon, zlat, zdvol, zdsal, zdtem, ingb)
	682	END DO
[11553]	683	END DO
	684	!
[11931]	685	!==============================================================================
	686	! 4.0: finalisation and compute ssh equivalent of the volume correction
	687	!==============================================================================
	688	!
[11553]	689	! mask (>0 out)
	690	risfcpl_cons_vol(:,:,: ) = risfcpl_cons_vol(:,:,: ) * tmask(:,:,:)
	691	risfcpl_cons_tsc(:,:,:,jp_sal) = risfcpl_cons_tsc(:,:,:,jp_sal) * tmask(:,:,:)
	692	risfcpl_cons_tsc(:,:,:,jp_tem) = risfcpl_cons_tsc(:,:,:,jp_tem) * tmask(:,:,:)
	693	!
	694	! add lbclnk
[14644]	695	CALL lbc_lnk( 'iscplrst', risfcpl_cons_tsc(:,:,:,jp_tem), 'T', 1.0_wp, risfcpl_cons_tsc(:,:,:,jp_sal), 'T', 1.0_wp, &
	696	& risfcpl_cons_vol(:,:,:) , 'T', 1.0_wp)
[11553]	697	!
	698	! ssh correction (for dynspg_ts)
	699	DO jk = 1,jpk
	700	risfcpl_cons_ssh(:,:) = risfcpl_cons_ssh(:,:) + risfcpl_cons_vol(:,:,jk)
	701	END DO
	702	risfcpl_cons_ssh(:,:) = risfcpl_cons_ssh(:,:) * r1_e1e2t(:,:)
	703	!
	704	END SUBROUTINE isfcpl_cons
	705	!
	706	SUBROUTINE update_isfpts(sisfpts, kpts, ki, kj, kk, pdvol, pdsal, pdtem, pratio, kfind)
	707	!!---------------------------------------------------------------------
	708	!! * ROUTINE update_isfpts *
	709	!!
	710	!! ** Purpose : if a cell become dry, we need to put the corrective increment elsewhere
	711	!!
	712	!! ** Action : update the list of point
	713	!!
	714	!!----------------------------------------------------------------------
	715	!!----------------------------------------------------------------------
	716	TYPE(isfcons), DIMENSION(:), INTENT(inout) :: sisfpts
[11931]	717	INTEGER, INTENT(inout) :: kpts
[11553]	718	!!----------------------------------------------------------------------
[14072]	719	INTEGER, INTENT(in ) :: ki, kj, kk ! target location (kfind=0)
[11931]	720	! ! or source location (kfind=1)
	721	INTEGER, INTENT(in ), OPTIONAL :: kfind ! 0 target cell already found
	722	! ! 1 target to be determined
[14072]	723	REAL(wp), INTENT(in ) :: pdvol, pdsal, pdtem, pratio ! vol/sal/tem increment
[11931]	724	! ! and ratio in case increment span over multiple cells.
[11553]	725	!!----------------------------------------------------------------------
	726	INTEGER :: ifind
	727	!!----------------------------------------------------------------------
[14072]	728	!
[11553]	729	! increment position
	730	kpts = kpts + 1
	731	!
	732	! define if we need to look for closest valid wet cell (no neighbours or neigbourg on halo)
	733	IF ( PRESENT(kfind) ) THEN
	734	ifind = kfind
	735	ELSE
	736	ifind = ( 1 - tmask_h(ki,kj) ) * tmask(ki,kj,kk)
	737	END IF
	738	!
	739	! update isfpts structure
	740	sisfpts(kpts) = isfcons(mig(ki), mjg(kj), kk, pratio * pdvol, pratio * pdsal, pratio * pdtem, glamt(ki,kj), gphit(ki,kj), ifind )
	741	!
	742	END SUBROUTINE update_isfpts
	743	!
[11931]	744	SUBROUTINE get_correction( ki, kj, kk, plon, plat, pvolinc, psalinc, pteminc, kfind)
[11553]	745	!!---------------------------------------------------------------------
	746	!! * ROUTINE get_correction *
	747	!!
	748	!! ** Action : - Find the closest valid cell if needed (wet and not on the halo)
	749	!! - Scale the correction depending of pratio (case where multiple wet neigbourgs)
	750	!! - Fill the correction array
	751	!!
	752	!!----------------------------------------------------------------------
[11931]	753	INTEGER , INTENT(in) :: ki, kj, kk, kfind ! target point indices
	754	REAL(wp), INTENT(in) :: plon, plat ! target point lon/lat
	755	REAL(wp), INTENT(in) :: pvolinc, pteminc,psalinc ! correction increment for vol/temp/salt
[11553]	756	!!----------------------------------------------------------------------
	757	INTEGER :: jj, ji, iig, ijg
	758	!!----------------------------------------------------------------------
	759	!
	760	! define global indice of correction location
	761	iig = ki ; ijg = kj
[11931]	762	IF ( kfind == 1 ) CALL dom_ngb( plon, plat, iig, ijg,'T', kk)
[11553]	763	!
	764	! fill the correction array
	765	DO jj = mj0(ijg),mj1(ijg)
	766	DO ji = mi0(iig),mi1(iig)
[11931]	767	! correct the vol_flx and corresponding heat/salt flx in the closest cell
[11553]	768	risfcpl_cons_vol(ji,jj,kk) = risfcpl_cons_vol(ji,jj,kk ) + pvolinc
	769	risfcpl_cons_tsc(ji,jj,kk,jp_sal) = risfcpl_cons_tsc(ji,jj,kk,jp_sal) + psalinc
	770	risfcpl_cons_tsc(ji,jj,kk,jp_tem) = risfcpl_cons_tsc(ji,jj,kk,jp_tem) + pteminc
	771	END DO
	772	END DO
	773
	774	END SUBROUTINE get_correction
	775
	776	END MODULE isfcpl

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

Download in other formats: