source: branches/UKMO/dev_r5518_GO6_package_inc_asm/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 8037

Last change on this file since 8037 was 8037, checked in by jwhile, 3 years ago

Further merge of CICE_with_keyasminc

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