New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcice_cice.F90 in branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 9124

Last change on this file since 9124 was 9124, checked in by gm, 6 years ago

dev_merge_2017: ln_timing instead of nn_timing + restricted timing to nemo_init and routine called by step in OPA_SRC

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