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

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

source: branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice_with_key_asminc/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 6175

Last change on this file since 6175 was 6175, checked in by kingr, 8 years ago

Possible fix for wroung dimension variables passed between NEMO and CICE.

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