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 NEMO/branches/UKMO/r12083_coupling_sequence/src/OCE/SBC – NEMO

source: NEMO/branches/UKMO/r12083_coupling_sequence/src/OCE/SBC/sbcice_cice.F90 @ 12465

Last change on this file since 12465 was 12465, checked in by jcastill, 2 years ago

Changes as in the original branch but updated to vn4.1

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