source: NEMO/branches/2020/ticket2487/src/OCE/SBC/sbcice_cice.F90 @ 13179

Last change on this file since 13179 was 13179, checked in by smueller, 3 months ago

Addition of initial sea-level compensation for non-zero sea-ice/snow mass (see ticket #2487)

  • Property svn:keywords set to Id
File size: 44.2 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 scalars
158      REAL(wp) ::   z1_area, zsshadj                  !   "      "
159      INTEGER  ::   ji, jj, jl, jk                    ! dummy loop indices
160      !!---------------------------------------------------------------------
161      !
162      IF(lwp) WRITE(numout,*)'cice_sbc_init'
163
164      ji_off = INT ( (jpiglo - nx_global) / 2 )
165      jj_off = INT ( (jpjglo - ny_global) / 2 )
166
167#if defined key_nemocice_decomp
168      ! Pass initial SST from NEMO to CICE so ice is initialised correctly if
169      ! there is no restart file.
170      ! Values from a CICE restart file would overwrite this
171      IF ( .NOT. ln_rstart ) THEN   
172         CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.) 
173      ENDIF 
174#endif
175
176! Initialize CICE
177      CALL CICE_Initialize
178
179! Do some CICE consistency checks
180      IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
181         IF ( calc_strair .OR. calc_Tsfc ) THEN
182            CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=F and calc_Tsfc=F in ice_in' )
183         ENDIF
184      ELSEIF (ksbc == jp_blk) THEN
185         IF ( .NOT. (calc_strair .AND. calc_Tsfc) ) THEN
186            CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=T and calc_Tsfc=T in ice_in' )
187         ENDIF
188      ENDIF
189
190
191! allocate sbc_ice and sbc_cice arrays
192      IF( sbc_ice_alloc()      /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate arrays' )
193      IF( sbc_ice_cice_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate cice arrays' )
194
195! Ensure ocean temperatures are nowhere below freezing if not a NEMO restart
196      IF( .NOT. ln_rstart ) THEN
197         tsn(:,:,:,jp_tem) = MAX (tsn(:,:,:,jp_tem),Tocnfrz)
198         tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
199      ENDIF
200
201      fr_iu(:,:)=0.0
202      fr_iv(:,:)=0.0
203
204      CALL cice2nemo(aice,fr_i, 'T', 1. )
205      IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
206         DO jl=1,ncat
207            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
208         ENDDO
209      ENDIF
210
211! T point to U point
212! T point to V point
213      DO jj=1,jpjm1
214         DO ji=1,jpim1
215            fr_iu(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji+1,jj))*umask(ji,jj,1)
216            fr_iv(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji,jj+1))*vmask(ji,jj,1)
217         ENDDO
218      ENDDO
219
220      CALL lbc_lnk_multi( 'sbcice_cice', fr_iu , 'U', 1.,  fr_iv , 'V', 1. )
221
222      ! set the snow+ice mass
223      CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
224      CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
225      snwice_mass  (:,:) = ( rhos * ztmp1(:,:) + rhoi * ztmp2(:,:)  )
226      snwice_mass_b(:,:) = snwice_mass(:,:)
227
228      IF( .NOT.ln_rstart ) THEN
229         IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
230            sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
231            sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
232
233!!gm This should be put elsewhere....   (same remark for limsbc)
234!!gm especially here it is assumed zstar coordinate, but it can be ztilde....
235            IF( .NOT.ln_linssh ) THEN
236               !
237               DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
238                  e3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
239                  e3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
240               ENDDO
241               e3t_a(:,:,:) = e3t_b(:,:,:)
242               ! Reconstruction of all vertical scale factors at now and before time-steps
243               ! =============================================================================
244               ! Horizontal scale factor interpolations
245               ! --------------------------------------
246               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
247               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
248               CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
249               CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
250               CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
251               ! Vertical scale factor interpolations
252               ! ------------------------------------
253               CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
254               CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
255               CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
256               CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
257               CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
258               ! t- and w- points depth
259               ! ----------------------
260               gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
261               gdepw_n(:,:,1) = 0.0_wp
262               gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
263               DO jk = 2, jpk
264                  gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk)
265                  gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
266                  gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn   (:,:)
267               END DO
268            ENDIF
269         ELSE
270            z1_area = 1.0_wp / glob_sum( 'sbcice_cice', e1e2t(:,:) )
271            zsshadj = glob_sum( 'sbcice_cice', e1e2t(:,:) * snwice_mass(:,:) ) * r1_rau0 * z1_area
272            WRITE(ctmp1,'(A36,F10.6,A24)') 'sbcice_cice:   mean ssh adjusted by ', -1.0_wp * zsshadj, ' m to compensate for the'
273            CALL ctl_warn( ctmp1,          '               initial ice+snow mass' )
274            sshn(:,:) = sshn(:,:) - zsshadj
275            sshb(:,:) = sshb(:,:) - zsshadj
276         ENDIF
277      ENDIF
278      !
279   END SUBROUTINE cice_sbc_init
280
281   
282   SUBROUTINE cice_sbc_in( kt, ksbc )
283      !!---------------------------------------------------------------------
284      !!                    ***  ROUTINE cice_sbc_in  ***
285      !! ** Purpose: Set coupling fields and pass to CICE
286      !!---------------------------------------------------------------------
287      INTEGER, INTENT(in   ) ::   kt   ! ocean time step
288      INTEGER, INTENT(in   ) ::   ksbc ! surface forcing type
289      !
290      INTEGER  ::   ji, jj, jl                   ! dummy loop indices     
291      REAL(wp), DIMENSION(jpi,jpj) :: ztmp, zpice
292      REAL(wp), DIMENSION(jpi,jpj,ncat) :: ztmpn
293      REAL(wp) ::   zintb, zintn  ! dummy argument
294      !!---------------------------------------------------------------------
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(:,:) ) * rLsub
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   END SUBROUTINE cice_sbc_in
496
497
498   SUBROUTINE cice_sbc_out( kt, ksbc )
499      !!---------------------------------------------------------------------
500      !!                    ***  ROUTINE cice_sbc_out  ***
501      !! ** Purpose: Get fields from CICE and set surface fields for NEMO
502      !!---------------------------------------------------------------------
503      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
504      INTEGER, INTENT( in  ) ::   ksbc ! surface forcing type
505     
506      INTEGER  ::   ji, jj, jl                 ! dummy loop indices
507      REAL(wp), DIMENSION(jpi,jpj) :: ztmp1, ztmp2
508      !!---------------------------------------------------------------------
509      !
510      IF( kt == nit000 )  THEN
511         IF(lwp) WRITE(numout,*)'cice_sbc_out'
512      ENDIF
513     
514! x comp of ocean-ice stress
515      CALL cice2nemo(strocnx,ztmp1,'F', -1. )
516      ss_iou(:,:)=0.0
517! F point to U point
518      DO jj=2,jpjm1
519         DO ji=2,jpim1
520            ss_iou(ji,jj) = 0.5 * ( ztmp1(ji,jj-1) + ztmp1(ji,jj) ) * umask(ji,jj,1)
521         ENDDO
522      ENDDO
523      CALL lbc_lnk( 'sbcice_cice', ss_iou , 'U', -1. )
524
525! y comp of ocean-ice stress
526      CALL cice2nemo(strocny,ztmp1,'F', -1. )
527      ss_iov(:,:)=0.0
528! F point to V point
529
530      DO jj=1,jpjm1
531         DO ji=2,jpim1
532            ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1)
533         ENDDO
534      ENDDO
535      CALL lbc_lnk( 'sbcice_cice', ss_iov , 'V', -1. )
536
537! x and y comps of surface stress
538! Combine wind stress and ocean-ice stress
539! [Note that fr_iu hasn't yet been updated, so still from start of CICE timestep]
540! strocnx and strocny already weighted by ice fraction in CICE so not done here
541
542      utau(:,:)=(1.0-fr_iu(:,:))*utau(:,:)-ss_iou(:,:)
543      vtau(:,:)=(1.0-fr_iv(:,:))*vtau(:,:)-ss_iov(:,:)     
544 
545! Also need ice/ocean stress on T points so that taum can be updated
546! This interpolation is already done in CICE so best to use those values
547      CALL cice2nemo(strocnxT,ztmp1,'T',-1.) 
548      CALL cice2nemo(strocnyT,ztmp2,'T',-1.) 
549 
550! Update taum with modulus of ice-ocean stress
551! strocnxT and strocnyT are not weighted by ice fraction in CICE so must be done here
552taum(:,:)=(1.0-fr_i(:,:))*taum(:,:)+fr_i(:,:)*SQRT(ztmp1*ztmp1 + ztmp2*ztmp2) 
553
554! Freshwater fluxes
555
556      IF (ksbc == jp_flx) THEN
557! Note that emp from the forcing files is evap*(1-aice)-(tprecip-aice*sprecip)
558! What we want here is evap*(1-aice)-tprecip*(1-aice) hence manipulation below
559! Not ideal since aice won't be the same as in the atmosphere. 
560! Better to use evap and tprecip? (but for now don't read in evap in this case)
561         emp(:,:)  = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
562      ELSE IF (ksbc == jp_blk) THEN
563         emp(:,:)  = (1.0-fr_i(:,:))*emp(:,:)       
564      ELSE IF (ksbc == jp_purecpl) THEN
565! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above)
566! This is currently as required with the coupling fields from the UM atmosphere
567         emp(:,:) = emp_tot(:,:)+tprecip(:,:)*fr_i(:,:) 
568      ENDIF
569
570#if defined key_cice4
571      CALL cice2nemo(fresh_gbm,ztmp1,'T', 1. )
572      CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1. )
573#else
574      CALL cice2nemo(fresh_ai,ztmp1,'T', 1. )
575      CALL cice2nemo(fsalt_ai,ztmp2,'T', 1. )
576#endif
577
578! Check to avoid unphysical expression when ice is forming (ztmp1 negative)
579! Otherwise we are effectively allowing ice of higher salinity than the ocean to form
580! which has to be compensated for by the ocean salinity potentially going negative
581! This check breaks conservation but seems reasonable until we have prognostic ice salinity
582! Note the 1000.0 below is to convert from kg salt to g salt (needed for PSU)
583      WHERE (ztmp1(:,:).lt.0.0) ztmp2(:,:)=MAX(ztmp2(:,:),ztmp1(:,:)*sss_m(:,:)/1000.0)
584      sfx(:,:)=ztmp2(:,:)*1000.0
585      emp(:,:)=emp(:,:)-ztmp1(:,:)
586      fmmflx(:,:) = ztmp1(:,:) !!Joakim edit
587     
588      CALL lbc_lnk_multi( 'sbcice_cice', emp , 'T', 1., sfx , 'T', 1. )
589
590! Solar penetrative radiation and non solar surface heat flux
591
592! Scale qsr and qns according to ice fraction (bulk formulae only)
593
594      IF (ksbc == jp_blk) THEN
595         qsr(:,:)=qsr(:,:)*(1.0-fr_i(:,:))
596         qns(:,:)=qns(:,:)*(1.0-fr_i(:,:))
597      ENDIF
598! Take into account snow melting except for fully coupled when already in qns_tot
599      IF (ksbc == jp_purecpl) THEN
600         qsr(:,:)= qsr_tot(:,:)
601         qns(:,:)= qns_tot(:,:)
602      ELSE
603         qns(:,:)= qns(:,:)-sprecip(:,:)*Lfresh*(1.0-fr_i(:,:))
604      ENDIF
605
606! Now add in ice / snow related terms
607! [fswthru will be zero unless running with calc_Tsfc=T in CICE]
608#if defined key_cice4
609      CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1. )
610#else
611      CALL cice2nemo(fswthru_ai,ztmp1,'T', 1. )
612#endif
613      qsr(:,:)=qsr(:,:)+ztmp1(:,:)
614      CALL lbc_lnk( 'sbcice_cice', qsr , 'T', 1. )
615
616      DO jj=1,jpj
617         DO ji=1,jpi
618            nfrzmlt(ji,jj)=MAX(nfrzmlt(ji,jj),0.0)
619         ENDDO
620      ENDDO
621
622#if defined key_cice4
623      CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1. )
624#else
625      CALL cice2nemo(fhocn_ai,ztmp1,'T', 1. )
626#endif
627      qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
628
629      CALL lbc_lnk( 'sbcice_cice', qns , 'T', 1. )
630
631! Prepare for the following CICE time-step
632
633      CALL cice2nemo(aice,fr_i,'T', 1. )
634      IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN
635         DO jl=1,ncat
636            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
637         ENDDO
638      ENDIF
639
640! T point to U point
641! T point to V point
642      DO jj=1,jpjm1
643         DO ji=1,jpim1
644            fr_iu(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji+1,jj))*umask(ji,jj,1)
645            fr_iv(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji,jj+1))*vmask(ji,jj,1)
646         ENDDO
647      ENDDO
648
649      CALL lbc_lnk_multi( 'sbcice_cice', fr_iu , 'U', 1., fr_iv , 'V', 1. )
650
651      ! set the snow+ice mass
652      CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
653      CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
654      snwice_mass  (:,:) = ( rhos * ztmp1(:,:) + rhoi * ztmp2(:,:)  )
655      snwice_mass_b(:,:) = snwice_mass(:,:)
656      snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt
657      !
658   END SUBROUTINE cice_sbc_out
659
660
661   SUBROUTINE cice_sbc_hadgam( kt )
662      !!---------------------------------------------------------------------
663      !!                    ***  ROUTINE cice_sbc_hadgam  ***
664      !! ** Purpose: Prepare fields needed to pass to HadGAM3 atmosphere
665      !!
666      !!
667      !!---------------------------------------------------------------------
668      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
669      !!
670      INTEGER  ::   jl                        ! dummy loop index
671      INTEGER  ::   ierror
672      !!---------------------------------------------------------------------
673      !
674      IF( kt == nit000 )  THEN
675         IF(lwp) WRITE(numout,*)'cice_sbc_hadgam'
676         IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
677      ENDIF
678
679      !                                         ! =========================== !
680      !                                         !   Prepare Coupling fields   !
681      !                                         ! =========================== !
682      !
683      ! x and y comp of ice velocity
684      !
685      CALL cice2nemo(uvel,u_ice,'F', -1. )
686      CALL cice2nemo(vvel,v_ice,'F', -1. )
687      !
688      ! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out 
689      !
690      ! Snow and ice thicknesses (CO_2 and CO_3)
691      !
692      DO jl = 1, ncat
693         CALL cice2nemo( vsnon(:,:,jl,:), h_s(:,:,jl),'T', 1. )
694         CALL cice2nemo( vicen(:,:,jl,:), h_i(:,:,jl),'T', 1. )
695      END DO
696      !
697   END SUBROUTINE cice_sbc_hadgam
698
699
700   SUBROUTINE cice_sbc_final
701      !!---------------------------------------------------------------------
702      !!                    ***  ROUTINE cice_sbc_final  ***
703      !! ** Purpose: Finalize CICE
704      !!---------------------------------------------------------------------
705      !
706      IF(lwp) WRITE(numout,*)'cice_sbc_final'
707      !
708      CALL CICE_Finalize
709      !
710   END SUBROUTINE cice_sbc_final
711
712
713   SUBROUTINE cice_sbc_force (kt)
714      !!---------------------------------------------------------------------
715      !!                    ***  ROUTINE cice_sbc_force  ***
716      !! ** Purpose : Provide CICE forcing from files
717      !!
718      !!---------------------------------------------------------------------
719      !! ** Method  :   READ monthly flux file in NetCDF files
720      !!     
721      !!  snowfall   
722      !!  rainfall   
723      !!  sublimation rate   
724      !!  topmelt (category)
725      !!  botmelt (category)
726      !!
727      !! History :
728      !!----------------------------------------------------------------------
729      USE iom
730      !!
731      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
732      !!
733      INTEGER  ::   ierror             ! return error code
734      INTEGER  ::   ifpr               ! dummy loop index
735      !!
736      CHARACTER(len=100) ::  cn_dir                            !   Root directory for location of CICE forcing files
737      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read
738      TYPE(FLD_N) ::   sn_snow, sn_rain, sn_sblm               ! informations about the fields to be read
739      TYPE(FLD_N) ::   sn_top1, sn_top2, sn_top3, sn_top4, sn_top5
740      TYPE(FLD_N) ::   sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5 
741      !!
742      NAMELIST/namsbc_cice/ cn_dir, sn_snow, sn_rain, sn_sblm,   &
743         &                          sn_top1, sn_top2, sn_top3, sn_top4, sn_top5,  &
744         &                          sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
745      INTEGER :: ios
746      !!---------------------------------------------------------------------
747
748      !                                         ! ====================== !
749      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
750         !                                      ! ====================== !
751         ! namsbc_cice is not yet in the reference namelist
752         ! set file information (default values)
753         cn_dir = './'       ! directory in which the model is executed
754
755         ! (NB: frequency positive => hours, negative => months)
756         !            !    file          ! frequency !  variable    ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   ! landmask
757         !            !    name          !  (hours)  !   name       !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      ! file
758         sn_snow = FLD_N( 'snowfall_1m'  ,    -1.    ,  'snowfall'  ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    ) 
759         sn_rain = FLD_N( 'rainfall_1m'  ,    -1.    ,  'rainfall'  ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    ) 
760         sn_sblm = FLD_N( 'sublim_1m'    ,    -1.    ,  'sublim'    ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
761         sn_top1 = FLD_N( 'topmeltn1_1m' ,    -1.    ,  'topmeltn1' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
762         sn_top2 = FLD_N( 'topmeltn2_1m' ,    -1.    ,  'topmeltn2' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
763         sn_top3 = FLD_N( 'topmeltn3_1m' ,    -1.    ,  'topmeltn3' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
764         sn_top4 = FLD_N( 'topmeltn4_1m' ,    -1.    ,  'topmeltn4' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
765         sn_top5 = FLD_N( 'topmeltn5_1m' ,    -1.    ,  'topmeltn5' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
766         sn_bot1 = FLD_N( 'botmeltn1_1m' ,    -1.    ,  'botmeltn1' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
767         sn_bot2 = FLD_N( 'botmeltn2_1m' ,    -1.    ,  'botmeltn2' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
768         sn_bot3 = FLD_N( 'botmeltn3_1m' ,    -1.    ,  'botmeltn3' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
769         sn_bot4 = FLD_N( 'botmeltn4_1m' ,    -1.    ,  'botmeltn4' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
770         sn_bot5 = FLD_N( 'botmeltn5_1m' ,    -1.    ,  'botmeltn5' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ,  ''    )
771
772         REWIND( numnam_ref )              ! Namelist namsbc_cice in reference namelist :
773         READ  ( numnam_ref, namsbc_cice, IOSTAT = ios, ERR = 901)
774901      IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_cice in reference namelist' )
775
776         REWIND( numnam_cfg )              ! Namelist namsbc_cice in configuration namelist : Parameters of the run
777         READ  ( numnam_cfg, namsbc_cice, IOSTAT = ios, ERR = 902 )
778902      IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_cice in configuration namelist' )
779         IF(lwm) WRITE ( numond, namsbc_cice )
780
781         ! store namelist information in an array
782         slf_i(jp_snow) = sn_snow   ;   slf_i(jp_rain) = sn_rain   ;   slf_i(jp_sblm) = sn_sblm
783         slf_i(jp_top1) = sn_top1   ;   slf_i(jp_top2) = sn_top2   ;   slf_i(jp_top3) = sn_top3
784         slf_i(jp_top4) = sn_top4   ;   slf_i(jp_top5) = sn_top5   ;   slf_i(jp_bot1) = sn_bot1
785         slf_i(jp_bot2) = sn_bot2   ;   slf_i(jp_bot3) = sn_bot3   ;   slf_i(jp_bot4) = sn_bot4
786         slf_i(jp_bot5) = sn_bot5
787         
788         ! set sf structure
789         ALLOCATE( sf(jpfld), STAT=ierror )
790         IF( ierror > 0 ) THEN
791            CALL ctl_stop( 'cice_sbc_force: unable to allocate sf structure' )   ;   RETURN
792         ENDIF
793
794         DO ifpr= 1, jpfld
795            ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
796            ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
797         END DO
798
799         ! fill sf with slf_i and control print
800         CALL fld_fill( sf, slf_i, cn_dir, 'cice_sbc_force', 'flux formulation for CICE', 'namsbc_cice' )
801         !
802      ENDIF
803
804      CALL fld_read( kt, nn_fsbc, sf )           ! Read input fields and provides the
805      !                                          ! input fields at the current time-step
806
807      ! set the fluxes from read fields
808      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1)
809      tprecip(:,:) = sf(jp_snow)%fnow(:,:,1)+sf(jp_rain)%fnow(:,:,1)
810! May be better to do this conversion somewhere else
811      qla_ice(:,:,1) = -rLsub*sf(jp_sblm)%fnow(:,:,1)
812      topmelt(:,:,1) = sf(jp_top1)%fnow(:,:,1)
813      topmelt(:,:,2) = sf(jp_top2)%fnow(:,:,1)
814      topmelt(:,:,3) = sf(jp_top3)%fnow(:,:,1)
815      topmelt(:,:,4) = sf(jp_top4)%fnow(:,:,1)
816      topmelt(:,:,5) = sf(jp_top5)%fnow(:,:,1)
817      botmelt(:,:,1) = sf(jp_bot1)%fnow(:,:,1)
818      botmelt(:,:,2) = sf(jp_bot2)%fnow(:,:,1)
819      botmelt(:,:,3) = sf(jp_bot3)%fnow(:,:,1)
820      botmelt(:,:,4) = sf(jp_bot4)%fnow(:,:,1)
821      botmelt(:,:,5) = sf(jp_bot5)%fnow(:,:,1)
822
823      ! control print (if less than 100 time-step asked)
824      IF( nitend-nit000 <= 100 .AND. lwp ) THEN
825         WRITE(numout,*) 
826         WRITE(numout,*) '        read forcing fluxes for CICE OK'
827         CALL FLUSH(numout)
828      ENDIF
829
830   END SUBROUTINE cice_sbc_force
831
832   SUBROUTINE nemo2cice( pn, pc, cd_type, psgn)
833      !!---------------------------------------------------------------------
834      !!                    ***  ROUTINE nemo2cice  ***
835      !! ** Purpose :   Transfer field in NEMO array to field in CICE array. 
836#if defined key_nemocice_decomp
837      !!             
838      !!                NEMO and CICE PE sub domains are identical, hence
839      !!                there is no need to gather or scatter data from
840      !!                one PE configuration to another.
841#else
842      !!                Automatically gather/scatter between
843      !!                different processors and blocks
844      !! ** Method :    A. Ensure all haloes are filled in NEMO field (pn)
845      !!                B. Gather pn into global array (png)
846      !!                C. Map png into CICE global array (pcg)
847      !!                D. Scatter pcg to CICE blocks (pc) + update haloes 
848#endif
849      !!---------------------------------------------------------------------
850      CHARACTER(len=1), INTENT( in ) ::   &
851          cd_type       ! nature of pn grid-point
852          !             !   = T or F gridpoints
853      REAL(wp), INTENT( in ) ::   &
854          psgn          ! control of the sign change
855          !             !   =-1 , the sign is modified following the type of b.c. used
856          !             !   = 1 , no sign change
857      REAL(wp), DIMENSION(jpi,jpj) :: pn
858#if !defined key_nemocice_decomp
859      REAL(wp), DIMENSION(jpiglo,jpjglo) :: png2
860      REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
861#endif
862      REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
863      INTEGER (int_kind) :: &
864         field_type,        &! id for type of field (scalar, vector, angle)
865         grid_loc            ! id for location on horizontal grid
866                            !  (center, NEcorner, Nface, Eface)
867
868      INTEGER  ::   ji, jj, jn                      ! dummy loop indices
869      !!---------------------------------------------------------------------
870
871!     A. Ensure all haloes are filled in NEMO field (pn)
872
873      CALL lbc_lnk( 'sbcice_cice', pn , cd_type, psgn )
874
875#if defined key_nemocice_decomp
876
877      ! Copy local domain data from NEMO to CICE field
878      pc(:,:,1)=0.0
879      DO jj=2,ny_block-1
880         DO ji=2,nx_block-1
881            pc(ji,jj,1)=pn(ji-1+ji_off,jj-1+jj_off)
882         ENDDO
883      ENDDO
884
885#else
886
887!     B. Gather pn into global array (png)
888
889      IF ( jpnij > 1) THEN
890         CALL mppsync
891         CALL mppgather (pn,0,png) 
892         CALL mppsync
893      ELSE
894         png(:,:,1)=pn(:,:)
895      ENDIF
896
897!     C. Map png into CICE global array (pcg)
898
899! Need to make sure this is robust to changes in NEMO halo rows....
900! (may be OK but not 100% sure)
901
902      IF (nproc==0) THEN     
903!        pcg(:,:)=0.0
904         DO jn=1,jpnij
905            DO jj=nldjt(jn),nlejt(jn)
906               DO ji=nldit(jn),nleit(jn)
907                  png2(ji+nimppt(jn)-1,jj+njmppt(jn)-1)=png(ji,jj,jn)
908               ENDDO
909            ENDDO
910         ENDDO
911         DO jj=1,ny_global
912            DO ji=1,nx_global
913               pcg(ji,jj)=png2(ji+ji_off,jj+jj_off)
914            ENDDO
915         ENDDO
916      ENDIF
917
918#endif
919
920      SELECT CASE ( cd_type )
921         CASE ( 'T' )
922            grid_loc=field_loc_center
923         CASE ( 'F' )                             
924            grid_loc=field_loc_NEcorner
925      END SELECT
926
927      SELECT CASE ( NINT(psgn) )
928         CASE ( -1 )
929            field_type=field_type_vector
930         CASE ( 1 )                             
931            field_type=field_type_scalar
932      END SELECT
933
934#if defined key_nemocice_decomp
935      ! Ensure CICE halos are up to date
936      CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
937#else
938!     D. Scatter pcg to CICE blocks (pc) + update halos
939      CALL scatter_global(pc, pcg, 0, distrb_info, grid_loc, field_type)
940#endif
941
942   END SUBROUTINE nemo2cice
943
944   SUBROUTINE cice2nemo ( pc, pn, cd_type, psgn )
945      !!---------------------------------------------------------------------
946      !!                    ***  ROUTINE cice2nemo  ***
947      !! ** Purpose :   Transfer field in CICE array to field in NEMO array.
948#if defined key_nemocice_decomp
949      !!             
950      !!                NEMO and CICE PE sub domains are identical, hence
951      !!                there is no need to gather or scatter data from
952      !!                one PE configuration to another.
953#else 
954      !!                Automatically deal with scatter/gather between
955      !!                different processors and blocks
956      !! ** Method :    A. Gather CICE blocks (pc) into global array (pcg)
957      !!                B. Map pcg into NEMO global array (png)
958      !!                C. Scatter png into NEMO field (pn) for each processor
959      !!                D. Ensure all haloes are filled in pn
960#endif
961      !!---------------------------------------------------------------------
962
963      CHARACTER(len=1), INTENT( in ) ::   &
964          cd_type       ! nature of pn grid-point
965          !             !   = T or F gridpoints
966      REAL(wp), INTENT( in ) ::   &
967          psgn          ! control of the sign change
968          !             !   =-1 , the sign is modified following the type of b.c. used
969          !             !   = 1 , no sign change
970      REAL(wp), DIMENSION(jpi,jpj) :: pn
971
972#if defined key_nemocice_decomp
973      INTEGER (int_kind) :: &
974         field_type,        & ! id for type of field (scalar, vector, angle)
975         grid_loc             ! id for location on horizontal grid
976                              ! (center, NEcorner, Nface, Eface)
977#else
978      REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
979#endif
980
981      REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
982
983      INTEGER  ::   ji, jj, jn                      ! dummy loop indices
984
985
986#if defined key_nemocice_decomp
987
988      SELECT CASE ( cd_type )
989         CASE ( 'T' )
990            grid_loc=field_loc_center
991         CASE ( 'F' )                             
992            grid_loc=field_loc_NEcorner
993      END SELECT
994
995      SELECT CASE ( NINT(psgn) )
996         CASE ( -1 )
997            field_type=field_type_vector
998         CASE ( 1 )                             
999            field_type=field_type_scalar
1000      END SELECT
1001
1002      CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
1003
1004
1005      pn(:,:)=0.0
1006      DO jj=1,jpjm1
1007         DO ji=1,jpim1
1008            pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1)
1009         ENDDO
1010      ENDDO
1011
1012#else
1013
1014!      A. Gather CICE blocks (pc) into global array (pcg)
1015
1016      CALL gather_global(pcg, pc, 0, distrb_info)
1017
1018!     B. Map pcg into NEMO global array (png)
1019
1020! Need to make sure this is robust to changes in NEMO halo rows....
1021! (may be OK but not spent much time thinking about it)
1022! Note that non-existent pcg elements may be used below, but
1023! the lbclnk call on pn will replace these with sensible values
1024
1025      IF (nproc==0) THEN
1026         png(:,:,:)=0.0
1027         DO jn=1,jpnij
1028            DO jj=nldjt(jn),nlejt(jn)
1029               DO ji=nldit(jn),nleit(jn)
1030                  png(ji,jj,jn)=pcg(ji+nimppt(jn)-1-ji_off,jj+njmppt(jn)-1-jj_off)
1031               ENDDO
1032            ENDDO
1033         ENDDO
1034      ENDIF
1035
1036!     C. Scatter png into NEMO field (pn) for each processor
1037
1038      IF ( jpnij > 1) THEN
1039         CALL mppsync
1040         CALL mppscatter (png,0,pn) 
1041         CALL mppsync
1042      ELSE
1043         pn(:,:)=png(:,:,1)
1044      ENDIF
1045
1046#endif
1047
1048!     D. Ensure all haloes are filled in pn
1049
1050      CALL lbc_lnk( 'sbcice_cice', pn , cd_type, psgn )
1051
1052   END SUBROUTINE cice2nemo
1053
1054#else
1055   !!----------------------------------------------------------------------
1056   !!   Default option           Dummy module         NO CICE sea-ice model
1057   !!----------------------------------------------------------------------
1058CONTAINS
1059
1060   SUBROUTINE sbc_ice_cice ( kt, ksbc )     ! Dummy routine
1061      IMPLICIT NONE
1062      INTEGER, INTENT( in ) :: kt, ksbc
1063      WRITE(*,*) 'sbc_ice_cice: You should not have seen this print! error?', kt
1064   END SUBROUTINE sbc_ice_cice
1065
1066   SUBROUTINE cice_sbc_init (ksbc)    ! Dummy routine
1067      IMPLICIT NONE
1068      INTEGER, INTENT( in ) :: ksbc
1069      WRITE(*,*) 'cice_sbc_init: You should not have seen this print! error?', ksbc
1070   END SUBROUTINE cice_sbc_init
1071
1072   SUBROUTINE cice_sbc_final     ! Dummy routine
1073      IMPLICIT NONE
1074      WRITE(*,*) 'cice_sbc_final: You should not have seen this print! error?'
1075   END SUBROUTINE cice_sbc_final
1076
1077#endif
1078
1079   !!======================================================================
1080END MODULE sbcice_cice
Note: See TracBrowser for help on using the repository browser.