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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_r8356_plus_form_drag/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 8484

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

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

Download in other formats: