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

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

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 8324

Last change on this file since 8324 was 8324, checked in by clem, 7 years ago

STEP3 (2): clean separation between sea-ice and ocean

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