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/SWE – NEMO

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/SWE/sbcice_cice.F90 @ 12983

Last change on this file since 12983 was 12983, checked in by techene, 4 years ago

SWE r12822 from Antoine updated with qco and rk3 options : draft 1

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