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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcice_cice.F90 in NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC – NEMO

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/SBC/sbcice_cice.F90 @ 12583

Last change on this file since 12583 was 12583, checked in by techene, 2 years ago

OCE/DOM/domqe.F90: add gdep at time level Kbb in dom_qe_sf_update, OCE/DOM/domzgr_substitute.h90: create the substitute module, OCE/DYN/dynatfLF.F90, OCE/TRA/traatfLF.F90: move boundary condition management and agrif management from atf modules to OCE/steplf.F90, OCE/SBC/sbcice_cice.F90, ICE/iceistate.F90 : remove dom_vvl_interpol and replace by dom_vvl_zgr ?

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