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/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 9819

Last change on this file since 9819 was 9819, checked in by dancopsey, 6 years ago
Added a bunch of fixes and a nemo2cice print statement.
File size: 46.9 KB

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

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

Download in other formats: