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/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 5987

Last change on this file since 5987 was 5987, checked in by timgraham, 8 years ago

Merged head of trunk (r5936) into branch

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