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

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

source: branches/UKMO/dev_r5518_nemo2cice_prints/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 9818

Last change on this file since 9818 was 9818, checked in by dancopsey, 6 years ago

Add min and max prints to nemo2cice.

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