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

Context Navigation

source: branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 4880

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

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

Download in other formats: