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/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5026

Last change on this file since 5026 was 5026, checked in by timgraham, 9 years ago

First set of changes for multicategory coupling with CICE

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