Ticket #1796: sbcice_cice.F90

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