source: branches/UKMO/dev_r5518_FOAMv14_output_heat_fluxes/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 10900

Last change on this file since 10900 was 10900, checked in by anaguiar, 2 years ago

Changes to allow checking qns balance and if all terms are correctly evaluated - Test 4 described in https://code.metoffice.gov.uk/trac/gmed/ticket/454#comment:2

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