source: branches/UKMO/dev_r5518_GO6_package_r8356_plus_form_drag/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 8819

Last change on this file since 8819 was 8819, checked in by jamrae, 3 years ago

Implemented code to pass sea ice skin roughness length to atmosphere.

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