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

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

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/SBC/sbcice_cice.F90 @ 13630

Last change on this file since 13630 was 13630, checked in by mocavero, 4 years ago

Add neighborhood collectives calls in the NEMO src - ticket #2496

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