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

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

source: branches/UKMO/dev_merge_2017_GC_couple_pkg/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 9723

Last change on this file since 9723 was 9723, checked in by dancopsey, 6 years ago

Removed call top cice_sbc_hadgam is this is being done in the dev_merge_2017_CICE_interface branch. This means that we must always use the dev_merge_2017_CICE_interface branch with this branch.

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