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.

sbcice_cice.F90 in branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 3625

Last change on this file since 3625 was 3625, checked in by acc, 12 years ago
Branch dev_NOC_2012_r3555. #1006. Step 7. Check in code now merged with dev_r3385_NOCS04_HAMF
File size: 41.9 KB

Rev	Line
[2874]	1	MODULE sbcice_cice
	2	!!======================================================================
	3	!! * MODULE sbcice_cice *
	4	!! To couple with sea ice model CICE (LANL)
	5	!!=====================================================================
	6	#if defined key_cice
	7	!!----------------------------------------------------------------------
	8	!! 'key_cice' : CICE sea-ice model
	9	!!----------------------------------------------------------------------
	10	!! sbc_ice_cice : sea-ice model time-stepping and update ocean sbc over ice-covered area
	11	!!
	12	!!
	13	!!----------------------------------------------------------------------
	14	USE oce ! ocean dynamics and tracers
	15	USE dom_oce ! ocean space and time domain
[3275]	16	USE domvvl
[3625]	17	USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic
[2874]	18	USE in_out_manager ! I/O manager
	19	USE lib_mpp ! distributed memory computing library
	20	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[3186]	21	USE wrk_nemo ! work arrays
[3193]	22	USE timing ! Timing
[2874]	23	USE daymod ! calendar
	24	USE fldread ! read input fields
	25
	26	USE sbc_oce ! Surface boundary condition: ocean fields
	27	USE sbc_ice ! Surface boundary condition: ice fields
	28	USE sbcblk_core ! Surface boundary condition: CORE bulk
	29	USE sbccpl
	30
	31	USE ice_kinds_mod
	32	USE ice_blocks
	33	USE ice_domain
	34	USE ice_domain_size
	35	USE ice_boundary
	36	USE ice_constants
	37	USE ice_gather_scatter
	38	USE ice_calendar, only: dt
[3625]	39	USE ice_state, only: aice,aicen,uvel,vvel,vsno,vsnon,vice,vicen
[2874]	40	USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow, &
[3189]	41	sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_gbm, &
	42	fresh_gbm,fhocn_gbm,fswthru_gbm,frzmlt, &
[2874]	43	flatn_f,fsurfn_f,fcondtopn_f, &
	44	uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl, &
	45	swvdr,swvdf,swidr,swidf
	46	USE ice_forcing, only: frcvdr,frcvdf,frcidr,frcidf
[3176]	47	USE ice_atmo, only: calc_strair
	48	USE ice_therm_vertical, only: calc_Tsfc
[2874]	49
	50	USE CICE_InitMod
	51	USE CICE_RunMod
	52	USE CICE_FinalMod
	53
	54	IMPLICIT NONE
	55	PRIVATE
	56
	57	!! * Routine accessibility
	58	PUBLIC cice_sbc_init ! routine called by sbc_init
	59	PUBLIC cice_sbc_final ! routine called by sbc_final
	60	PUBLIC sbc_ice_cice ! routine called by sbc
	61
[3625]	62	INTEGER , PARAMETER :: ji_off = INT ( (jpiglo - nx_global) / 2 )
	63	INTEGER , PARAMETER :: jj_off = INT ( (jpjglo - ny_global) / 2 )
	64
[2874]	65	INTEGER , PARAMETER :: jpfld = 13 ! maximum number of files to read
	66	INTEGER , PARAMETER :: jp_snow = 1 ! index of snow file
	67	INTEGER , PARAMETER :: jp_rain = 2 ! index of rain file
	68	INTEGER , PARAMETER :: jp_sblm = 3 ! index of sublimation file
	69	INTEGER , PARAMETER :: jp_top1 = 4 ! index of category 1 topmelt file
	70	INTEGER , PARAMETER :: jp_top2 = 5 ! index of category 2 topmelt file
	71	INTEGER , PARAMETER :: jp_top3 = 6 ! index of category 3 topmelt file
	72	INTEGER , PARAMETER :: jp_top4 = 7 ! index of category 4 topmelt file
	73	INTEGER , PARAMETER :: jp_top5 = 8 ! index of category 5 topmelt file
	74	INTEGER , PARAMETER :: jp_bot1 = 9 ! index of category 1 botmelt file
	75	INTEGER , PARAMETER :: jp_bot2 = 10 ! index of category 2 botmelt file
	76	INTEGER , PARAMETER :: jp_bot3 = 11 ! index of category 3 botmelt file
	77	INTEGER , PARAMETER :: jp_bot4 = 12 ! index of category 4 botmelt file
	78	INTEGER , PARAMETER :: jp_bot5 = 13 ! index of category 5 botmelt file
	79	TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read)
	80
	81	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:), PRIVATE :: png ! local array used in sbc_cice_ice
	82
	83	!! * Substitutions
	84	# include "domzgr_substitute.h90"
	85
	86	CONTAINS
	87
	88	INTEGER FUNCTION sbc_ice_cice_alloc()
	89	!!----------------------------------------------------------------------
	90	!! * FUNCTION sbc_ice_cice_alloc *
	91	!!----------------------------------------------------------------------
	92	ALLOCATE( png(jpi,jpj,jpnij), STAT=sbc_ice_cice_alloc )
	93	IF( lk_mpp ) CALL mpp_sum ( sbc_ice_cice_alloc )
	94	IF( sbc_ice_cice_alloc > 0 ) CALL ctl_warn('sbc_ice_cice_alloc: allocation of arrays failed.')
	95	END FUNCTION sbc_ice_cice_alloc
	96
	97	SUBROUTINE sbc_ice_cice( kt, nsbc )
	98	!!---------------------------------------------------------------------
	99	!! * ROUTINE sbc_ice_cice *
	100	!!
	101	!! ** Purpose : update the ocean surface boundary condition via the
	102	!! CICE Sea Ice Model time stepping
	103	!!
[3040]	104	!! ** Method : - Get any extra forcing fields for CICE
	105	!! - Prepare forcing fields
[2874]	106	!! - CICE model time stepping
	107	!! - call the routine that computes mass and
	108	!! heat fluxes at the ice/ocean interface
	109	!!
	110	!! ** Action : - time evolution of the CICE sea-ice model
	111	!! - update all sbc variables below sea-ice:
[3625]	112	!! utau, vtau, qns , qsr, emp , sfx
[2874]	113	!!---------------------------------------------------------------------
	114	INTEGER, INTENT(in) :: kt ! ocean time step
	115	INTEGER, INTENT(in) :: nsbc ! surface forcing type
	116	!!----------------------------------------------------------------------
[3193]	117	!
	118	IF( nn_timing == 1 ) CALL timing_start('sbc_ice_cice')
	119	!
[2874]	120	! !----------------------!
	121	IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN ! Ice time-step only !
	122	! !----------------------!
	123
	124	! Make sure any fluxes required for CICE are set
	125	IF ( nsbc == 2 ) THEN
	126	CALL cice_sbc_force(kt)
	127	ELSE IF ( nsbc == 5 ) THEN
	128	CALL sbc_cpl_ice_flx( 1.0-fr_i )
	129	ENDIF
	130
	131	CALL cice_sbc_in ( kt, nsbc )
	132	CALL CICE_Run
	133	CALL cice_sbc_out ( kt, nsbc )
	134
	135	IF ( nsbc == 5 ) CALL cice_sbc_hadgam(kt+1)
	136
	137	ENDIF ! End sea-ice time step only
[3193]	138	!
	139	IF( nn_timing == 1 ) CALL timing_stop('sbc_ice_cice')
[2874]	140
	141	END SUBROUTINE sbc_ice_cice
	142
	143	SUBROUTINE cice_sbc_init (nsbc)
	144	!!---------------------------------------------------------------------
	145	!! * ROUTINE cice_sbc_init *
[3040]	146	!! ** Purpose: Initialise ice related fields for NEMO and coupling
[2874]	147	!!
[3625]	148	INTEGER, INTENT( in ) :: nsbc ! surface forcing type
	149	REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
	150	REAL(wp) :: zcoefu, zcoefv, zcoeff ! local scalar
	151	INTEGER :: ji, jj, jl ! dummy loop indices
[2874]	152	!!---------------------------------------------------------------------
	153
[3193]	154	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_init')
	155	!
[3625]	156	CALL wrk_alloc( jpi,jpj, ztmp1, ztmp2 )
	157	!
[2874]	158	IF(lwp) WRITE(numout,*)'cice_sbc_init'
	159
	160	! Initialize CICE
[3176]	161	CALL CICE_Initialize
[2874]	162
[3176]	163	! Do some CICE consistency checks
[3193]	164	IF ( (nsbc == 2) .OR. (nsbc == 5) ) THEN
	165	IF ( calc_strair .OR. calc_Tsfc ) THEN
	166	CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=F and calc_Tsfc=F in ice_in' )
	167	ENDIF
	168	ELSEIF (nsbc == 4) THEN
	169	IF ( .NOT. (calc_strair .AND. calc_Tsfc) ) THEN
	170	CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=T and calc_Tsfc=T in ice_in' )
	171	ENDIF
	172	ENDIF
[3176]	173
	174
[2874]	175	! allocate sbc_ice and sbc_cice arrays
	176	IF( sbc_ice_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate arrays' )
	177	IF( sbc_ice_cice_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate cice arrays' )
	178
	179	! Ensure ocean temperatures are nowhere below freezing if not a NEMO restart
	180	IF( .NOT. ln_rstart ) THEN
	181	tsn(:,:,:,jp_tem) = MAX (tsn(:,:,:,jp_tem),Tocnfrz)
	182	tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
	183	ENDIF
	184
[3193]	185	fr_iu(:,:)=0.0
	186	fr_iv(:,:)=0.0
[2874]	187
[3193]	188	CALL cice2nemo(aice,fr_i, 'T', 1. )
	189	IF ( (nsbc == 2).OR.(nsbc == 5) ) THEN
[3625]	190	DO jl=1,ncat
	191	CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
[3193]	192	ENDDO
	193	ENDIF
[2874]	194
	195	! T point to U point
	196	! T point to V point
[3193]	197	DO jj=1,jpjm1
	198	DO ji=1,jpim1
	199	fr_iu(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji+1,jj))umask(ji,jj,1)
	200	fr_iv(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji,jj+1))vmask(ji,jj,1)
	201	ENDDO
	202	ENDDO
[2874]	203
[3193]	204	CALL lbc_lnk ( fr_iu , 'U', 1. )
	205	CALL lbc_lnk ( fr_iv , 'V', 1. )
[3625]	206
	207	! ! embedded sea ice
	208	IF( nn_ice_embd /= 0 ) THEN ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
	209	CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
	210	CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
	211	snwice_mass (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:) )
	212	snwice_mass_b(:,:) = snwice_mass(:,:)
	213	ELSE
	214	snwice_mass (:,:) = 0.0_wp ! no mass exchanges
	215	snwice_mass_b(:,:) = 0.0_wp ! no mass exchanges
	216	ENDIF
	217	IF( nn_ice_embd == 2 .AND. & ! full embedment (case 2) & no restart :
	218	& .NOT.ln_rstart ) THEN ! deplete the initial ssh belew sea-ice area
	219	sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
	220	sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
	221	!
	222	! Note: Changed the initial values of sshb and sshn=> need to recompute ssh[u,v,f]_[b,n]
	223	! which were previously set in domvvl
	224	IF ( lk_vvl ) THEN ! Is this necessary? embd 2 should be restricted to vvl only???
	225	DO jj = 1, jpjm1
	226	DO ji = 1, jpim1 ! caution: use of Vector Opt. not possible
	227	zcoefu = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
	228	zcoefv = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
	229	zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
	230	sshu_b(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) &
	231	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
	232	sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) &
	233	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
	234	sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) &
	235	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
	236	sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) &
	237	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
	238	END DO
	239	END DO
	240	CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. )
	241	CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. )
	242	DO jj = 1, jpjm1
	243	DO ji = 1, jpim1 ! NO Vector Opt.
	244	sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) &
	245	& / ( e1f(ji,jj ) * e2f(ji,jj ) ) &
	246	& * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) &
	247	& + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
	248	END DO
	249	END DO
	250	CALL lbc_lnk( sshf_n, 'F', 1. )
	251	ENDIF
	252	ENDIF
	253
	254	CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
[3193]	255	!
	256	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_init')
	257	!
[2874]	258	END SUBROUTINE cice_sbc_init
	259
[3152]	260
[2874]	261	SUBROUTINE cice_sbc_in (kt, nsbc)
	262	!!---------------------------------------------------------------------
	263	!! * ROUTINE cice_sbc_in *
[3040]	264	!! ** Purpose: Set coupling fields and pass to CICE
[2874]	265	!!---------------------------------------------------------------------
[3152]	266	INTEGER, INTENT(in ) :: kt ! ocean time step
	267	INTEGER, INTENT(in ) :: nsbc ! surface forcing type
[2874]	268
[3625]	269	INTEGER :: ji, jj, jl ! dummy loop indices
	270	REAL(wp), DIMENSION(:,:), POINTER :: ztmp, zpice
[3152]	271	REAL(wp), DIMENSION(:,:,:), POINTER :: ztmpn
[3625]	272	REAL(wp) :: zintb, zintn ! dummy argument
[3152]	273	!!---------------------------------------------------------------------
[2874]	274
[3193]	275	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_in')
	276	!
[3625]	277	CALL wrk_alloc( jpi,jpj, ztmp, zpice )
[3152]	278	CALL wrk_alloc( jpi,jpj,ncat, ztmpn )
[2874]	279
[3193]	280	IF( kt == nit000 ) THEN
[2874]	281	IF(lwp) WRITE(numout,*)'cice_sbc_in'
[3193]	282	ENDIF
[2874]	283
[3193]	284	ztmp(:,:)=0.0
[2874]	285
	286	! Aggregate ice concentration already set in cice_sbc_out (or cice_sbc_init on
	287	! the first time-step)
	288
	289	! forced and coupled case
	290
[3193]	291	IF ( (nsbc == 2).OR.(nsbc == 5) ) THEN
[2874]	292
[3193]	293	ztmpn(:,:,:)=0.0
[2874]	294
	295	! x comp of wind stress (CI_1)
	296	! U point to F point
[3193]	297	DO jj=1,jpjm1
	298	DO ji=1,jpi
	299	ztmp(ji,jj) = 0.5 * ( fr_iu(ji,jj) * utau(ji,jj) &
	300	+ fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1)
	301	ENDDO
	302	ENDDO
	303	CALL nemo2cice(ztmp,strax,'F', -1. )
[2874]	304
	305	! y comp of wind stress (CI_2)
	306	! V point to F point
[3193]	307	DO jj=1,jpj
	308	DO ji=1,jpim1
	309	ztmp(ji,jj) = 0.5 * ( fr_iv(ji,jj) * vtau(ji,jj) &
	310	+ fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
	311	ENDDO
	312	ENDDO
	313	CALL nemo2cice(ztmp,stray,'F', -1. )
[2874]	314
	315	! Surface downward latent heat flux (CI_5)
[3193]	316	IF (nsbc == 2) THEN
[3625]	317	DO jl=1,ncat
	318	ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl)
[3193]	319	ENDDO
	320	ELSE
[2874]	321	! emp_ice is set in sbc_cpl_ice_flx as sublimation-snow
[3193]	322	qla_ice(:,:,1)= - ( emp_ice(:,:)+sprecip(:,:) ) * Lsub
[2874]	323	! End of temporary code
[3193]	324	DO jj=1,jpj
	325	DO ji=1,jpi
	326	IF (fr_i(ji,jj).eq.0.0) THEN
[3625]	327	DO jl=1,ncat
	328	ztmpn(ji,jj,jl)=0.0
[3193]	329	ENDDO
	330	! This will then be conserved in CICE
	331	ztmpn(ji,jj,1)=qla_ice(ji,jj,1)
	332	ELSE
[3625]	333	DO jl=1,ncat
	334	ztmpn(ji,jj,jl)=qla_ice(ji,jj,1)*a_i(ji,jj,jl)/fr_i(ji,jj)
[3193]	335	ENDDO
	336	ENDIF
	337	ENDDO
	338	ENDDO
	339	ENDIF
[3625]	340	DO jl=1,ncat
	341	CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. )
[2874]	342
	343	! GBM conductive flux through ice (CI_6)
	344	! Convert to GBM
[3193]	345	IF (nsbc == 2) THEN
[3625]	346	ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl)
[3193]	347	ELSE
[3625]	348	ztmp(:,:) = botmelt(:,:,jl)
[3193]	349	ENDIF
[3625]	350	CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1. )
[2874]	351
	352	! GBM surface heat flux (CI_7)
	353	! Convert to GBM
[3193]	354	IF (nsbc == 2) THEN
[3625]	355	ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl)
[3193]	356	ELSE
[3625]	357	ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))
[3193]	358	ENDIF
[3625]	359	CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. )
[3193]	360	ENDDO
[2874]	361
[3193]	362	ELSE IF (nsbc == 4) THEN
[2874]	363
	364	! Pass CORE forcing fields to CICE (which will calculate heat fluxes etc itself)
	365	! x comp and y comp of atmosphere surface wind (CICE expects on T points)
[3193]	366	ztmp(:,:) = wndi_ice(:,:)
	367	CALL nemo2cice(ztmp,uatm,'T', -1. )
	368	ztmp(:,:) = wndj_ice(:,:)
	369	CALL nemo2cice(ztmp,vatm,'T', -1. )
	370	ztmp(:,:) = SQRT ( wndi_ice(:,:)2 + wndj_ice(:,:)2 )
	371	CALL nemo2cice(ztmp,wind,'T', 1. ) ! Wind speed (m/s)
	372	ztmp(:,:) = qsr_ice(:,:,1)
	373	CALL nemo2cice(ztmp,fsw,'T', 1. ) ! Incoming short-wave (W/m^2)
	374	ztmp(:,:) = qlw_ice(:,:,1)
	375	CALL nemo2cice(ztmp,flw,'T', 1. ) ! Incoming long-wave (W/m^2)
	376	ztmp(:,:) = tatm_ice(:,:)
	377	CALL nemo2cice(ztmp,Tair,'T', 1. ) ! Air temperature (K)
	378	CALL nemo2cice(ztmp,potT,'T', 1. ) ! Potential temp (K)
[2874]	379	! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows
[3193]	380	ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )
	381	! Constant (101000.) atm pressure assumed
	382	CALL nemo2cice(ztmp,rhoa,'T', 1. ) ! Air density (kg/m^3)
	383	ztmp(:,:) = qatm_ice(:,:)
	384	CALL nemo2cice(ztmp,Qa,'T', 1. ) ! Specific humidity (kg/kg)
	385	ztmp(:,:)=10.0
	386	CALL nemo2cice(ztmp,zlvl,'T', 1. ) ! Atmos level height (m)
[2874]	387
	388	! May want to check all values are physically realistic (as in CICE routine
	389	! prepare_forcing)?
	390
	391	! Divide shortwave into spectral bands (as in prepare_forcing)
[3193]	392	ztmp(:,:)=qsr_ice(:,:,1)*frcvdr ! visible direct
[2874]	393	CALL nemo2cice(ztmp,swvdr,'T', 1. )
[3193]	394	ztmp(:,:)=qsr_ice(:,:,1)*frcvdf ! visible diffuse
[2874]	395	CALL nemo2cice(ztmp,swvdf,'T', 1. )
[3193]	396	ztmp(:,:)=qsr_ice(:,:,1)*frcidr ! near IR direct
[2874]	397	CALL nemo2cice(ztmp,swidr,'T', 1. )
[3193]	398	ztmp(:,:)=qsr_ice(:,:,1)*frcidf ! near IR diffuse
[2874]	399	CALL nemo2cice(ztmp,swidf,'T', 1. )
	400
	401	ENDIF
	402
	403	! Snowfall
	404	! Ensure fsnow is positive (as in CICE routine prepare_forcing)
[3193]	405	ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0)
	406	CALL nemo2cice(ztmp,fsnow,'T', 1. )
[2874]	407
	408	! Rainfall
[3193]	409	ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
	410	CALL nemo2cice(ztmp,frain,'T', 1. )
[2874]	411
	412	! Freezing/melting potential
[3275]	413	! Calculated over NEMO leapfrog timestep (hence 2*dt)
[3277]	414	nfrzmlt(:,:)=rau0rcpfse3t_m(:,:,1)(Tocnfrz-sst_m(:,:))/(2.0dt)
[2874]	415
[3193]	416	ztmp(:,:) = nfrzmlt(:,:)
	417	CALL nemo2cice(ztmp,frzmlt,'T', 1. )
[2874]	418
	419	! SST and SSS
	420
[3193]	421	CALL nemo2cice(sst_m,sst,'T', 1. )
	422	CALL nemo2cice(sss_m,sss,'T', 1. )
[2874]	423
	424	! x comp and y comp of surface ocean current
	425	! U point to F point
[3193]	426	DO jj=1,jpjm1
	427	DO ji=1,jpi
	428	ztmp(ji,jj)=0.5(ssu_m(ji,jj)+ssu_m(ji,jj+1))fmask(ji,jj,1)
	429	ENDDO
	430	ENDDO
	431	CALL nemo2cice(ztmp,uocn,'F', -1. )
[2874]	432
	433	! V point to F point
[3193]	434	DO jj=1,jpj
	435	DO ji=1,jpim1
	436	ztmp(ji,jj)=0.5(ssv_m(ji,jj)+ssv_m(ji+1,jj))fmask(ji,jj,1)
	437	ENDDO
	438	ENDDO
	439	CALL nemo2cice(ztmp,vocn,'F', -1. )
[2874]	440
[3625]	441	IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==!
	442	!
	443	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
	444	! = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
	445	zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
	446	!
	447	! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
	448	! = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
	449	zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
	450	!
	451	zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
	452	!
	453	!
	454	ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==!
	455	zpice(:,:) = ssh_m(:,:)
	456	ENDIF
	457
[3189]	458	! x comp and y comp of sea surface slope (on F points)
	459	! T point to F point
[3193]	460	DO jj=1,jpjm1
	461	DO ji=1,jpim1
[3625]	462	ztmp(ji,jj)=0.5 * ( (zpice(ji+1,jj )-zpice(ji,jj ))/e1u(ji,jj ) &
	463	+ (zpice(ji+1,jj+1)-zpice(ji,jj+1))/e1u(ji,jj+1) ) &
[3193]	464	* fmask(ji,jj,1)
	465	ENDDO
	466	ENDDO
	467	CALL nemo2cice(ztmp,ss_tltx,'F', -1. )
[3189]	468
	469	! T point to F point
[3193]	470	DO jj=1,jpjm1
	471	DO ji=1,jpim1
[3625]	472	ztmp(ji,jj)=0.5 * ( (zpice(ji ,jj+1)-zpice(ji ,jj))/e2v(ji ,jj) &
	473	+ (zpice(ji+1,jj+1)-zpice(ji+1,jj))/e2v(ji+1,jj) ) &
[3193]	474	* fmask(ji,jj,1)
	475	ENDDO
	476	ENDDO
	477	CALL nemo2cice(ztmp,ss_tlty,'F', -1. )
[3189]	478
[3152]	479	CALL wrk_dealloc( jpi,jpj, ztmp )
	480	CALL wrk_dealloc( jpi,jpj,ncat, ztmpn )
[3193]	481	!
	482	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_in')
	483	!
[2874]	484	END SUBROUTINE cice_sbc_in
	485
[3152]	486
[2874]	487	SUBROUTINE cice_sbc_out (kt,nsbc)
	488	!!---------------------------------------------------------------------
	489	!! * ROUTINE cice_sbc_out *
[3040]	490	!! ** Purpose: Get fields from CICE and set surface fields for NEMO
[3152]	491	!!---------------------------------------------------------------------
[2874]	492	INTEGER, INTENT( in ) :: kt ! ocean time step
	493	INTEGER, INTENT( in ) :: nsbc ! surface forcing type
[3152]	494
[3625]	495	INTEGER :: ji, jj, jl ! dummy loop indices
	496	REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
[2874]	497	!!---------------------------------------------------------------------
	498
[3193]	499	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_out')
	500	!
[3625]	501	CALL wrk_alloc( jpi,jpj, ztmp1, ztmp2 )
[3152]	502
	503	IF( kt == nit000 ) THEN
[2874]	504	IF(lwp) WRITE(numout,*)'cice_sbc_out'
[3152]	505	ENDIF
	506
[2874]	507	! x comp of ocean-ice stress
[3625]	508	CALL cice2nemo(strocnx,ztmp1,'F', -1. )
[3193]	509	ss_iou(:,:)=0.0
[2874]	510	! F point to U point
[3193]	511	DO jj=2,jpjm1
	512	DO ji=2,jpim1
[3625]	513	ss_iou(ji,jj) = 0.5 * ( ztmp1(ji,jj-1) + ztmp1(ji,jj) ) * umask(ji,jj,1)
[3193]	514	ENDDO
	515	ENDDO
	516	CALL lbc_lnk( ss_iou , 'U', -1. )
[2874]	517
	518	! y comp of ocean-ice stress
[3625]	519	CALL cice2nemo(strocny,ztmp1,'F', -1. )
[3193]	520	ss_iov(:,:)=0.0
[2874]	521	! F point to V point
	522
[3193]	523	DO jj=1,jpjm1
	524	DO ji=2,jpim1
[3625]	525	ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1)
[3193]	526	ENDDO
	527	ENDDO
	528	CALL lbc_lnk( ss_iov , 'V', -1. )
[2874]	529
	530	! x and y comps of surface stress
	531	! Combine wind stress and ocean-ice stress
	532	! [Note that fr_iu hasn't yet been updated, so still from start of CICE timestep]
	533
[3193]	534	utau(:,:)=(1.0-fr_iu(:,:))*utau(:,:)-ss_iou(:,:)
	535	vtau(:,:)=(1.0-fr_iv(:,:))*vtau(:,:)-ss_iov(:,:)
[2874]	536
	537	! Freshwater fluxes
	538
[3193]	539	IF (nsbc == 2) THEN
[2874]	540	! Note that emp from the forcing files is evap(1-aice)-(tprecip-aicesprecip)
	541	! What we want here is evap(1-aice)-tprecip(1-aice) hence manipulation below
	542	! Not ideal since aice won't be the same as in the atmosphere.
	543	! Better to use evap and tprecip? (but for now don't read in evap in this case)
[3193]	544	emp(:,:) = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
	545	ELSE IF (nsbc == 4) THEN
	546	emp(:,:) = (1.0-fr_i(:,:))*emp(:,:)
	547	ELSE IF (nsbc ==5) THEN
[3625]	548	! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above)
	549	! This is currently as required with the coupling fields from the UM atmosphere
[3193]	550	emp(:,:) = emp_tot(:,:)+tprecip(:,:)*fr_i(:,:)
	551	ENDIF
[2874]	552
[3625]	553	CALL cice2nemo(fresh_gbm,ztmp1,'T', 1. )
	554	CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1. )
[2874]	555
[3625]	556	! Check to avoid unphysical expression when ice is forming (ztmp1 negative)
	557	! Otherwise we are effectively allowing ice of higher salinity than the ocean to form
	558	! which has to be compensated for by the ocean salinity potentially going negative
	559	! This check breaks conservation but seems reasonable until we have prognostic ice salinity
	560	! Note the 1000.0 below is to convert from kg salt to g salt (needed for PSU)
	561	WHERE (ztmp1(:,:).lt.0.0) ztmp2(:,:)=MAX(ztmp2(:,:),ztmp1(:,:)*sss_m(:,:)/1000.0)
	562	sfx(:,:)=ztmp2(:,:)*1000.0
	563	emp(:,:)=emp(:,:)-ztmp1(:,:)
	564
[3193]	565	CALL lbc_lnk( emp , 'T', 1. )
[3625]	566	CALL lbc_lnk( sfx , 'T', 1. )
[2874]	567
	568	! Solar penetrative radiation and non solar surface heat flux
	569
	570	! Scale qsr and qns according to ice fraction (bulk formulae only)
	571
[3193]	572	IF (nsbc == 4) THEN
	573	qsr(:,:)=qsr(:,:)*(1.0-fr_i(:,:))
	574	qns(:,:)=qns(:,:)*(1.0-fr_i(:,:))
	575	ENDIF
[2874]	576	! Take into account snow melting except for fully coupled when already in qns_tot
[3193]	577	IF (nsbc == 5) THEN
	578	qsr(:,:)= qsr_tot(:,:)
	579	qns(:,:)= qns_tot(:,:)
	580	ELSE
	581	qns(:,:)= qns(:,:)-sprecip(:,:)Lfresh(1.0-fr_i(:,:))
	582	ENDIF
[2874]	583
	584	! Now add in ice / snow related terms
	585	! [fswthru will be zero unless running with calc_Tsfc=T in CICE]
[3625]	586	CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1. )
	587	qsr(:,:)=qsr(:,:)+ztmp1(:,:)
[3193]	588	CALL lbc_lnk( qsr , 'T', 1. )
[2874]	589
[3193]	590	DO jj=1,jpj
	591	DO ji=1,jpi
[2874]	592	nfrzmlt(ji,jj)=MAX(nfrzmlt(ji,jj),0.0)
[3193]	593	ENDDO
	594	ENDDO
[2874]	595
[3625]	596	CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1. )
	597	qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
[2874]	598
[3193]	599	CALL lbc_lnk( qns , 'T', 1. )
[2874]	600
	601	! Prepare for the following CICE time-step
	602
[3193]	603	CALL cice2nemo(aice,fr_i,'T', 1. )
	604	IF ( (nsbc == 2).OR.(nsbc == 5) ) THEN
[3625]	605	DO jl=1,ncat
	606	CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
[3193]	607	ENDDO
	608	ENDIF
[2874]	609
	610	! T point to U point
	611	! T point to V point
[3193]	612	DO jj=1,jpjm1
	613	DO ji=1,jpim1
	614	fr_iu(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji+1,jj))umask(ji,jj,1)
	615	fr_iv(ji,jj)=0.5(fr_i(ji,jj)+fr_i(ji,jj+1))vmask(ji,jj,1)
	616	ENDDO
	617	ENDDO
[2874]	618
[3193]	619	CALL lbc_lnk ( fr_iu , 'U', 1. )
	620	CALL lbc_lnk ( fr_iv , 'V', 1. )
[2874]	621
[3625]	622	! ! embedded sea ice
	623	IF( nn_ice_embd /= 0 ) THEN ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
	624	CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
	625	CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
	626	snwice_mass (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:) )
	627	snwice_mass_b(:,:) = snwice_mass(:,:)
	628	snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt
	629	ENDIF
	630
[2874]	631	! Release work space
	632
[3625]	633	CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
[3193]	634	!
	635	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_out')
	636	!
[2874]	637	END SUBROUTINE cice_sbc_out
	638
[3152]	639
[2874]	640	#if defined key_oasis3 \|\| defined key_oasis4
	641	SUBROUTINE cice_sbc_hadgam( kt )
	642	!!---------------------------------------------------------------------
	643	!! * ROUTINE cice_sbc_hadgam *
[3040]	644	!! ** Purpose: Prepare fields needed to pass to HadGAM3 atmosphere
[2874]	645	!!
	646	!!
	647	INTEGER, INTENT( in ) :: kt ! ocean time step
	648	!!---------------------------------------------------------------------
	649
[3625]	650	INTEGER :: jl ! dummy loop index
[3193]	651	INTEGER :: ierror
[2874]	652
[3193]	653	IF( nn_timing == 1 ) CALL timing_start('cice_sbc_hadgam')
	654	!
	655	IF( kt == nit000 ) THEN
[2874]	656	IF(lwp) WRITE(numout,*)'cice_sbc_hadgam'
	657	IF( sbc_cpl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
[3193]	658	ENDIF
[2874]	659
	660	! ! =========================== !
	661	! ! Prepare Coupling fields !
	662	! ! =========================== !
	663
	664	! x and y comp of ice velocity
	665
[3193]	666	CALL cice2nemo(uvel,u_ice,'F', -1. )
	667	CALL cice2nemo(vvel,v_ice,'F', -1. )
[2874]	668
	669	! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out
	670
	671	! Snow and ice thicknesses (CO_2 and CO_3)
	672
[3625]	673	DO jl = 1,ncat
	674	CALL cice2nemo(vsnon(:,:,jl,:),ht_s(:,:,jl),'T', 1. )
	675	CALL cice2nemo(vicen(:,:,jl,:),ht_i(:,:,jl),'T', 1. )
[3193]	676	ENDDO
	677	!
	678	IF( nn_timing == 1 ) CALL timing_stop('cice_sbc_hadgam')
	679	!
[2874]	680	END SUBROUTINE cice_sbc_hadgam
	681
	682	#else
[2884]	683	SUBROUTINE cice_sbc_hadgam( kt ) ! Dummy routine
	684	INTEGER, INTENT( in ) :: kt ! ocean time step
[2874]	685	WRITE(,) 'cice_sbc_hadgam: You should not have seen this print! error?'
	686	END SUBROUTINE cice_sbc_hadgam
	687	#endif
	688
	689	SUBROUTINE cice_sbc_final
	690	!!---------------------------------------------------------------------
	691	!! * ROUTINE cice_sbc_final *
	692	!! ** Purpose: Finalize CICE
	693	!!---------------------------------------------------------------------
	694
	695	IF(lwp) WRITE(numout,*)'cice_sbc_final'
	696
[3193]	697	CALL CICE_Finalize
[2874]	698
	699	END SUBROUTINE cice_sbc_final
	700
	701	SUBROUTINE cice_sbc_force (kt)
	702	!!---------------------------------------------------------------------
	703	!! * ROUTINE cice_sbc_force *
	704	!! ** Purpose : Provide CICE forcing from files
	705	!!
	706	!!---------------------------------------------------------------------
	707	!! ** Method : READ monthly flux file in NetCDF files
	708	!!
	709	!! snowfall
	710	!! rainfall
	711	!! sublimation rate
	712	!! topmelt (category)
	713	!! botmelt (category)
	714	!!
	715	!! History :
	716	!!----------------------------------------------------------------------
	717	!! * Modules used
	718	USE iom
	719
	720	!! * arguments
	721	INTEGER, INTENT( in ) :: kt ! ocean time step
	722
	723	INTEGER :: ierror ! return error code
	724	INTEGER :: ifpr ! dummy loop index
	725	!!
	726	CHARACTER(len=100) :: cn_dir ! Root directory for location of CICE forcing files
	727	TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read
	728	TYPE(FLD_N) :: sn_snow, sn_rain, sn_sblm ! informations about the fields to be read
	729	TYPE(FLD_N) :: sn_top1, sn_top2, sn_top3, sn_top4, sn_top5
	730	TYPE(FLD_N) :: sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
	731
	732	!!
	733	NAMELIST/namsbc_cice/ cn_dir, sn_snow, sn_rain, sn_sblm, &
	734	& sn_top1, sn_top2, sn_top3, sn_top4, sn_top5, &
	735	& sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
	736	!!---------------------------------------------------------------------
	737
	738	! ! ====================== !
	739	IF( kt == nit000 ) THEN ! First call kt=nit000 !
	740	! ! ====================== !
	741	! set file information (default values)
	742	cn_dir = './' ! directory in which the model is executed
	743
	744	! (NB: frequency positive => hours, negative => months)
	745	! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation !
	746	! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs !
	747	sn_snow = FLD_N( 'snowfall_1m' , -1. , 'snowfall' , .true. , .true. , ' yearly' , '' , '' )
	748	sn_rain = FLD_N( 'rainfall_1m' , -1. , 'rainfall' , .true. , .true. , ' yearly' , '' , '' )
	749	sn_sblm = FLD_N( 'sublim_1m' , -1. , 'sublim' , .true. , .true. , ' yearly' , '' , '' )
	750	sn_top1 = FLD_N( 'topmeltn1_1m' , -1. , 'topmeltn1' , .true. , .true. , ' yearly' , '' , '' )
	751	sn_top2 = FLD_N( 'topmeltn2_1m' , -1. , 'topmeltn2' , .true. , .true. , ' yearly' , '' , '' )
	752	sn_top3 = FLD_N( 'topmeltn3_1m' , -1. , 'topmeltn3' , .true. , .true. , ' yearly' , '' , '' )
	753	sn_top4 = FLD_N( 'topmeltn4_1m' , -1. , 'topmeltn4' , .true. , .true. , ' yearly' , '' , '' )
	754	sn_top5 = FLD_N( 'topmeltn5_1m' , -1. , 'topmeltn5' , .true. , .true. , ' yearly' , '' , '' )
	755	sn_bot1 = FLD_N( 'botmeltn1_1m' , -1. , 'botmeltn1' , .true. , .true. , ' yearly' , '' , '' )
	756	sn_bot2 = FLD_N( 'botmeltn2_1m' , -1. , 'botmeltn2' , .true. , .true. , ' yearly' , '' , '' )
	757	sn_bot3 = FLD_N( 'botmeltn3_1m' , -1. , 'botmeltn3' , .true. , .true. , ' yearly' , '' , '' )
	758	sn_bot4 = FLD_N( 'botmeltn4_1m' , -1. , 'botmeltn4' , .true. , .true. , ' yearly' , '' , '' )
	759	sn_bot5 = FLD_N( 'botmeltn5_1m' , -1. , 'botmeltn5' , .true. , .true. , ' yearly' , '' , '' )
	760
	761	! REWIND ( numnam ) ! ... at some point might read in from NEMO namelist?
	762	! READ ( numnam, namsbc_cice )
	763
	764	! store namelist information in an array
	765	slf_i(jp_snow) = sn_snow ; slf_i(jp_rain) = sn_rain ; slf_i(jp_sblm) = sn_sblm
	766	slf_i(jp_top1) = sn_top1 ; slf_i(jp_top2) = sn_top2 ; slf_i(jp_top3) = sn_top3
	767	slf_i(jp_top4) = sn_top4 ; slf_i(jp_top5) = sn_top5 ; slf_i(jp_bot1) = sn_bot1
	768	slf_i(jp_bot2) = sn_bot2 ; slf_i(jp_bot3) = sn_bot3 ; slf_i(jp_bot4) = sn_bot4
	769	slf_i(jp_bot5) = sn_bot5
	770
	771	! set sf structure
	772	ALLOCATE( sf(jpfld), STAT=ierror )
	773	IF( ierror > 0 ) THEN
	774	CALL ctl_stop( 'cice_sbc_force: unable to allocate sf structure' ) ; RETURN
	775	ENDIF
	776
	777	DO ifpr= 1, jpfld
	778	ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
	779	ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
	780	END DO
	781
	782	! fill sf with slf_i and control print
	783	CALL fld_fill( sf, slf_i, cn_dir, 'cice_sbc_force', 'flux formulation for CICE', 'namsbc_cice' )
	784	!
	785	ENDIF
	786
	787	CALL fld_read( kt, nn_fsbc, sf ) ! Read input fields and provides the
	788	! ! input fields at the current time-step
	789
	790	! set the fluxes from read fields
	791	sprecip(:,:) = sf(jp_snow)%fnow(:,:,1)
	792	tprecip(:,:) = sf(jp_snow)%fnow(:,:,1)+sf(jp_rain)%fnow(:,:,1)
[3040]	793	! May be better to do this conversion somewhere else
[2874]	794	qla_ice(:,:,1) = -Lsub*sf(jp_sblm)%fnow(:,:,1)
	795	topmelt(:,:,1) = sf(jp_top1)%fnow(:,:,1)
	796	topmelt(:,:,2) = sf(jp_top2)%fnow(:,:,1)
	797	topmelt(:,:,3) = sf(jp_top3)%fnow(:,:,1)
	798	topmelt(:,:,4) = sf(jp_top4)%fnow(:,:,1)
	799	topmelt(:,:,5) = sf(jp_top5)%fnow(:,:,1)
	800	botmelt(:,:,1) = sf(jp_bot1)%fnow(:,:,1)
	801	botmelt(:,:,2) = sf(jp_bot2)%fnow(:,:,1)
	802	botmelt(:,:,3) = sf(jp_bot3)%fnow(:,:,1)
	803	botmelt(:,:,4) = sf(jp_bot4)%fnow(:,:,1)
	804	botmelt(:,:,5) = sf(jp_bot5)%fnow(:,:,1)
	805
	806	! control print (if less than 100 time-step asked)
	807	IF( nitend-nit000 <= 100 .AND. lwp ) THEN
	808	WRITE(numout,*)
	809	WRITE(numout,*) ' read forcing fluxes for CICE OK'
	810	CALL FLUSH(numout)
	811	ENDIF
	812
	813	END SUBROUTINE cice_sbc_force
	814
	815	SUBROUTINE nemo2cice( pn, pc, cd_type, psgn)
	816	!!---------------------------------------------------------------------
	817	!! * ROUTINE nemo2cice *
	818	!! ** Purpose : Transfer field in NEMO array to field in CICE array.
	819	#if defined key_nemocice_decomp
	820	!!
	821	!! NEMO and CICE PE sub domains are identical, hence
	822	!! there is no need to gather or scatter data from
	823	!! one PE configuration to another.
	824	#else
	825	!! Automatically gather/scatter between
	826	!! different processors and blocks
	827	!! ** Method : A. Ensure all haloes are filled in NEMO field (pn)
	828	!! B. Gather pn into global array (png)
	829	!! C. Map png into CICE global array (pcg)
	830	!! D. Scatter pcg to CICE blocks (pc) + update haloes
	831	#endif
	832	!!---------------------------------------------------------------------
	833
[3193]	834	CHARACTER(len=1), INTENT( in ) :: &
	835	cd_type ! nature of pn grid-point
	836	! ! = T or F gridpoints
	837	REAL(wp), INTENT( in ) :: &
	838	psgn ! control of the sign change
	839	! ! =-1 , the sign is modified following the type of b.c. used
	840	! ! = 1 , no sign change
	841	REAL(wp), DIMENSION(jpi,jpj) :: pn
[2874]	842	#if !defined key_nemocice_decomp
[3625]	843	REAL(wp), DIMENSION(jpiglo,jpjglo) :: png2
[3193]	844	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
[2874]	845	#endif
[3193]	846	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
	847	INTEGER (int_kind) :: &
	848	field_type, &! id for type of field (scalar, vector, angle)
	849	grid_loc ! id for location on horizontal grid
[2874]	850	! (center, NEcorner, Nface, Eface)
	851
[3193]	852	INTEGER :: ji, jj, jn ! dummy loop indices
[2874]	853
[3193]	854	! A. Ensure all haloes are filled in NEMO field (pn)
[2874]	855
[3193]	856	CALL lbc_lnk( pn , cd_type, psgn )
[2874]	857
	858	#if defined key_nemocice_decomp
	859
[3193]	860	! Copy local domain data from NEMO to CICE field
	861	pc(:,:,1)=0.0
[3625]	862	DO jj=2,ny_block-1
	863	DO ji=2,nx_block-1
	864	pc(ji,jj,1)=pn(ji-1+ji_off,jj-1+jj_off)
[3193]	865	ENDDO
	866	ENDDO
[2874]	867
	868	#else
	869
[3193]	870	! B. Gather pn into global array (png)
[2874]	871
[3193]	872	IF ( jpnij > 1) THEN
	873	CALL mppsync
	874	CALL mppgather (pn,0,png)
	875	CALL mppsync
	876	ELSE
	877	png(:,:,1)=pn(:,:)
	878	ENDIF
[2874]	879
[3193]	880	! C. Map png into CICE global array (pcg)
[2874]	881
	882	! Need to make sure this is robust to changes in NEMO halo rows....
	883	! (may be OK but not 100% sure)
	884
[3193]	885	IF (nproc==0) THEN
[2874]	886	! pcg(:,:)=0.0
[3193]	887	DO jn=1,jpnij
[3625]	888	DO jj=nldjt(jn),nlejt(jn)
	889	DO ji=nldit(jn),nleit(jn)
	890	png2(ji+nimppt(jn)-1,jj+njmppt(jn)-1)=png(ji,jj,jn)
[3193]	891	ENDDO
	892	ENDDO
	893	ENDDO
[3625]	894	DO jj=1,ny_global
	895	DO ji=1,nx_global
	896	pcg(ji,jj)=png2(ji+ji_off,jj+jj_off)
	897	ENDDO
	898	ENDDO
[3193]	899	ENDIF
[2874]	900
	901	#endif
	902
[3193]	903	SELECT CASE ( cd_type )
	904	CASE ( 'T' )
	905	grid_loc=field_loc_center
	906	CASE ( 'F' )
	907	grid_loc=field_loc_NEcorner
	908	END SELECT
[2874]	909
[3193]	910	SELECT CASE ( NINT(psgn) )
	911	CASE ( -1 )
	912	field_type=field_type_vector
	913	CASE ( 1 )
	914	field_type=field_type_scalar
	915	END SELECT
[2874]	916
	917	#if defined key_nemocice_decomp
[3193]	918	! Ensure CICE halos are up to date
	919	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
[2874]	920	#else
[3193]	921	! D. Scatter pcg to CICE blocks (pc) + update halos
	922	CALL scatter_global(pc, pcg, 0, distrb_info, grid_loc, field_type)
[2874]	923	#endif
	924
	925	END SUBROUTINE nemo2cice
	926
	927	SUBROUTINE cice2nemo ( pc, pn, cd_type, psgn )
	928	!!---------------------------------------------------------------------
	929	!! * ROUTINE cice2nemo *
	930	!! ** Purpose : Transfer field in CICE array to field in NEMO array.
	931	#if defined key_nemocice_decomp
	932	!!
	933	!! NEMO and CICE PE sub domains are identical, hence
	934	!! there is no need to gather or scatter data from
	935	!! one PE configuration to another.
	936	#else
	937	!! Automatically deal with scatter/gather between
	938	!! different processors and blocks
	939	!! ** Method : A. Gather CICE blocks (pc) into global array (pcg)
	940	!! B. Map pcg into NEMO global array (png)
	941	!! C. Scatter png into NEMO field (pn) for each processor
	942	!! D. Ensure all haloes are filled in pn
	943	#endif
	944	!!---------------------------------------------------------------------
	945
[3193]	946	CHARACTER(len=1), INTENT( in ) :: &
	947	cd_type ! nature of pn grid-point
	948	! ! = T or F gridpoints
	949	REAL(wp), INTENT( in ) :: &
	950	psgn ! control of the sign change
	951	! ! =-1 , the sign is modified following the type of b.c. used
	952	! ! = 1 , no sign change
	953	REAL(wp), DIMENSION(jpi,jpj) :: pn
[2874]	954
	955	#if defined key_nemocice_decomp
[3193]	956	INTEGER (int_kind) :: &
	957	field_type, & ! id for type of field (scalar, vector, angle)
	958	grid_loc ! id for location on horizontal grid
	959	! (center, NEcorner, Nface, Eface)
[2874]	960	#else
[3193]	961	REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
[2874]	962	#endif
	963
[3193]	964	REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
[2874]	965
[3193]	966	INTEGER :: ji, jj, jn ! dummy loop indices
[2874]	967
	968
	969	#if defined key_nemocice_decomp
	970
[3193]	971	SELECT CASE ( cd_type )
	972	CASE ( 'T' )
	973	grid_loc=field_loc_center
	974	CASE ( 'F' )
	975	grid_loc=field_loc_NEcorner
	976	END SELECT
[2874]	977
[3193]	978	SELECT CASE ( NINT(psgn) )
	979	CASE ( -1 )
	980	field_type=field_type_vector
	981	CASE ( 1 )
	982	field_type=field_type_scalar
	983	END SELECT
[2874]	984
[3193]	985	CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
[2874]	986
	987
[3193]	988	pn(:,:)=0.0
	989	DO jj=1,jpjm1
	990	DO ji=1,jpim1
[3625]	991	pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1)
[3193]	992	ENDDO
	993	ENDDO
[2874]	994
	995	#else
	996
[3193]	997	! A. Gather CICE blocks (pc) into global array (pcg)
[2874]	998
[3193]	999	CALL gather_global(pcg, pc, 0, distrb_info)
[2874]	1000
	1001	! B. Map pcg into NEMO global array (png)
	1002
	1003	! Need to make sure this is robust to changes in NEMO halo rows....
	1004	! (may be OK but not spent much time thinking about it)
[3625]	1005	! Note that non-existent pcg elements may be used below, but
	1006	! the lbclnk call on pn will replace these with sensible values
[2874]	1007
[3193]	1008	IF (nproc==0) THEN
	1009	png(:,:,:)=0.0
	1010	DO jn=1,jpnij
[3625]	1011	DO jj=nldjt(jn),nlejt(jn)
	1012	DO ji=nldit(jn),nleit(jn)
	1013	png(ji,jj,jn)=pcg(ji+nimppt(jn)-1-ji_off,jj+njmppt(jn)-1-jj_off)
[3193]	1014	ENDDO
	1015	ENDDO
	1016	ENDDO
	1017	ENDIF
[2874]	1018
[3193]	1019	! C. Scatter png into NEMO field (pn) for each processor
[2874]	1020
[3193]	1021	IF ( jpnij > 1) THEN
	1022	CALL mppsync
	1023	CALL mppscatter (png,0,pn)
	1024	CALL mppsync
	1025	ELSE
	1026	pn(:,:)=png(:,:,1)
	1027	ENDIF
[2874]	1028
	1029	#endif
	1030
[3193]	1031	! D. Ensure all haloes are filled in pn
[2874]	1032
[3193]	1033	CALL lbc_lnk( pn , cd_type, psgn )
[2874]	1034
	1035	END SUBROUTINE cice2nemo
	1036
	1037	#else
	1038	!!----------------------------------------------------------------------
	1039	!! Default option Dummy module NO CICE sea-ice model
	1040	!!----------------------------------------------------------------------
	1041	CONTAINS
	1042
	1043	SUBROUTINE sbc_ice_cice ( kt, nsbc ) ! Dummy routine
	1044	WRITE(,) 'sbc_ice_cice: You should not have seen this print! error?', kt
	1045	END SUBROUTINE sbc_ice_cice
	1046
	1047	SUBROUTINE cice_sbc_init (nsbc) ! Dummy routine
	1048	WRITE(,) 'cice_sbc_init: You should not have seen this print! error?'
	1049	END SUBROUTINE cice_sbc_init
	1050
	1051	SUBROUTINE cice_sbc_final ! Dummy routine
	1052	WRITE(,) 'cice_sbc_final: You should not have seen this print! error?'
	1053	END SUBROUTINE cice_sbc_final
	1054
	1055	#endif
	1056
	1057	!!======================================================================
	1058	END MODULE sbcice_cice

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

Download in other formats: