source: branches/UKMO/dev_merge_2017_CICE_interface/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 9695

Last change on this file since 9695 was 9695, checked in by dancopsey, 2 years ago

Apply -Lsub to all cases of zevap_ice. Multiply v_ip by a_i.

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