source: branches/UKMO/dev_r5518_GO6_package_STARTHOUR/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 11274

Last change on this file since 11274 was 11274, checked in by dancopsey, 14 months ago

Merge in sublimation bug fixes (GMED#449).

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