source: branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 10276

Last change on this file since 10276 was 10276, checked in by emmafiedler, 23 months ago

Freeboard assimilation updates

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