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/2012/dev_3352_UKMO8_CICE/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2012/dev_3352_UKMO8_CICE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 3355

Last change on this file since 3355 was 3355, checked in by charris, 12 years ago

#953 First set of modifications (details in ticket).

File size: 37.6 KB
Line 
1MODULE sbcice_cice
2   !!======================================================================
3   !!                       ***  MODULE  sbcice_cice  ***
4   !! To couple with sea ice model CICE (LANL)
5   !!=====================================================================
6#if defined key_cice
7   !!----------------------------------------------------------------------
8   !!   'key_cice' :                                     CICE sea-ice model
9   !!----------------------------------------------------------------------
10   !!   sbc_ice_cice  : sea-ice model time-stepping and update ocean sbc over ice-covered area
11   !!   
12   !!   
13   !!----------------------------------------------------------------------
14   USE oce             ! ocean dynamics and tracers
15   USE dom_oce         ! ocean space and time domain
16   USE domvvl
17   USE phycst, only : rcp, rau0
18   USE in_out_manager  ! I/O manager
19   USE lib_mpp         ! distributed memory computing library
20   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
21   USE wrk_nemo        ! work arrays
22   USE timing          ! Timing
23   USE daymod          ! calendar
24   USE fldread         ! read input fields
25
26   USE sbc_oce         ! Surface boundary condition: ocean fields
27   USE sbc_ice         ! Surface boundary condition: ice   fields
28   USE sbcblk_core     ! Surface boundary condition: CORE bulk
29   USE sbccpl
30
31   USE ice_kinds_mod
32   USE ice_blocks
33   USE ice_domain
34   USE ice_domain_size
35   USE ice_boundary
36   USE ice_constants
37   USE ice_gather_scatter
38   USE ice_calendar, only: dt
39   USE ice_state, only: aice,aicen,uvel,vvel,vsnon,vicen
40   USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow,  &
41                sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_gbm,     &
42                fresh_gbm,fhocn_gbm,fswthru_gbm,frzmlt,          &
43                flatn_f,fsurfn_f,fcondtopn_f,                    &
44                uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl,   &
45                swvdr,swvdf,swidr,swidf
46   USE ice_forcing, only: frcvdr,frcvdf,frcidr,frcidf
47   USE ice_atmo, only: calc_strair
48   USE ice_therm_vertical, only: calc_Tsfc
49
50   USE CICE_InitMod
51   USE CICE_RunMod
52   USE CICE_FinalMod
53
54   IMPLICIT NONE
55   PRIVATE
56
57   !! * Routine accessibility
58   PUBLIC cice_sbc_init   ! routine called by sbc_init
59   PUBLIC cice_sbc_final  ! routine called by sbc_final
60   PUBLIC sbc_ice_cice    ! routine called by sbc
61
62   INTEGER , PARAMETER ::   jpfld   = 13   ! maximum number of files to read
63   INTEGER , PARAMETER ::   jp_snow = 1    ! index of snow file
64   INTEGER , PARAMETER ::   jp_rain = 2    ! index of rain file
65   INTEGER , PARAMETER ::   jp_sblm = 3    ! index of sublimation file
66   INTEGER , PARAMETER ::   jp_top1 = 4    ! index of category 1 topmelt file
67   INTEGER , PARAMETER ::   jp_top2 = 5    ! index of category 2 topmelt file
68   INTEGER , PARAMETER ::   jp_top3 = 6    ! index of category 3 topmelt file
69   INTEGER , PARAMETER ::   jp_top4 = 7    ! index of category 4 topmelt file
70   INTEGER , PARAMETER ::   jp_top5 = 8    ! index of category 5 topmelt file
71   INTEGER , PARAMETER ::   jp_bot1 = 9    ! index of category 1 botmelt file
72   INTEGER , PARAMETER ::   jp_bot2 = 10   ! index of category 2 botmelt file
73   INTEGER , PARAMETER ::   jp_bot3 = 11   ! index of category 3 botmelt file
74   INTEGER , PARAMETER ::   jp_bot4 = 12   ! index of category 4 botmelt file
75   INTEGER , PARAMETER ::   jp_bot5 = 13   ! index of category 5 botmelt file
76   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read)
77
78   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:), PRIVATE ::   png     ! local array used in sbc_cice_ice
79
80   !! * Substitutions
81#  include "domzgr_substitute.h90"
82
83CONTAINS
84
85   INTEGER FUNCTION sbc_ice_cice_alloc()
86      !!----------------------------------------------------------------------
87      !!                ***  FUNCTION sbc_ice_cice_alloc  ***
88      !!----------------------------------------------------------------------
89      ALLOCATE( png(jpi,jpj,jpnij), STAT=sbc_ice_cice_alloc )
90      IF( lk_mpp                 )   CALL mpp_sum ( sbc_ice_cice_alloc )
91      IF( sbc_ice_cice_alloc > 0 )   CALL ctl_warn('sbc_ice_cice_alloc: allocation of arrays failed.')
92   END FUNCTION sbc_ice_cice_alloc
93
94   SUBROUTINE sbc_ice_cice( kt, nsbc )
95      !!---------------------------------------------------------------------
96      !!                  ***  ROUTINE sbc_ice_cice  ***
97      !!                   
98      !! ** Purpose :   update the ocean surface boundary condition via the
99      !!                CICE Sea Ice Model time stepping
100      !!
101      !! ** Method  : - Get any extra forcing fields for CICE 
102      !!              - Prepare forcing fields
103      !!              - CICE model time stepping
104      !!              - call the routine that computes mass and
105      !!                heat fluxes at the ice/ocean interface
106      !!
107      !! ** Action  : - time evolution of the CICE sea-ice model
108      !!              - update all sbc variables below sea-ice:
109      !!                utau, vtau, qns , qsr, emp , emps
110      !!---------------------------------------------------------------------
111      INTEGER, INTENT(in) ::   kt      ! ocean time step
112      INTEGER, INTENT(in) ::   nsbc    ! surface forcing type
113      !!----------------------------------------------------------------------
114      !
115      IF( nn_timing == 1 )  CALL timing_start('sbc_ice_cice')
116      !
117      !                                        !----------------------!
118      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  !
119         !                                     !----------------------!
120
121         ! Make sure any fluxes required for CICE are set
122         IF ( nsbc == 2 )  THEN
123            CALL cice_sbc_force(kt)
124         ELSE IF ( nsbc == 5 ) THEN
125            CALL sbc_cpl_ice_flx( 1.0-fr_i  )
126         ENDIF
127
128         CALL cice_sbc_in ( kt, nsbc )
129         CALL CICE_Run
130         CALL cice_sbc_out ( kt, nsbc )
131
132         IF ( nsbc == 5 )  CALL cice_sbc_hadgam(kt+1)
133
134      ENDIF                                          ! End sea-ice time step only
135      !
136      IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_cice')
137
138   END SUBROUTINE sbc_ice_cice
139
140   SUBROUTINE cice_sbc_init (nsbc)
141      !!---------------------------------------------------------------------
142      !!                    ***  ROUTINE cice_sbc_init  ***
143      !! ** Purpose: Initialise ice related fields for NEMO and coupling
144      !!
145      INTEGER, INTENT( in  ) ::   nsbc ! surface forcing type
146      !!---------------------------------------------------------------------
147
148      INTEGER  ::   ji, jj, jl                        ! dummy loop indices
149
150      IF( nn_timing == 1 )  CALL timing_start('cice_sbc_init')
151      !
152      IF(lwp) WRITE(numout,*)'cice_sbc_init'
153
154! Initialize CICE
155      CALL CICE_Initialize
156
157! Do some CICE consistency checks
158      IF ( (nsbc == 2) .OR. (nsbc == 5) ) THEN
159         IF ( calc_strair .OR. calc_Tsfc ) THEN
160            CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=F and calc_Tsfc=F in ice_in' )
161         ENDIF
162      ELSEIF (nsbc == 4) THEN
163         IF ( .NOT. (calc_strair .AND. calc_Tsfc) ) THEN
164            CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=T and calc_Tsfc=T in ice_in' )
165         ENDIF
166      ENDIF
167
168
169! allocate sbc_ice and sbc_cice arrays
170      IF( sbc_ice_alloc()      /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate arrays' )
171      IF( sbc_ice_cice_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate cice arrays' )
172
173! Ensure ocean temperatures are nowhere below freezing if not a NEMO restart
174      IF( .NOT. ln_rstart ) THEN
175         tsn(:,:,:,jp_tem) = MAX (tsn(:,:,:,jp_tem),Tocnfrz)
176         tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
177      ENDIF
178
179      fr_iu(:,:)=0.0
180      fr_iv(:,:)=0.0
181
182      CALL cice2nemo(aice,fr_i, 'T', 1. )
183      IF ( (nsbc == 2).OR.(nsbc == 5) ) THEN
184         DO jl=1,ncat
185            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
186         ENDDO
187      ENDIF
188
189! T point to U point
190! T point to V point
191      DO jj=1,jpjm1
192         DO ji=1,jpim1
193            fr_iu(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji+1,jj))*umask(ji,jj,1)
194            fr_iv(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji,jj+1))*vmask(ji,jj,1)
195         ENDDO
196      ENDDO
197
198      CALL lbc_lnk ( fr_iu , 'U', 1. )
199      CALL lbc_lnk ( fr_iv , 'V', 1. )
200      !
201      IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_init')
202      !
203   END SUBROUTINE cice_sbc_init
204
205   
206   SUBROUTINE cice_sbc_in (kt, nsbc)
207      !!---------------------------------------------------------------------
208      !!                    ***  ROUTINE cice_sbc_in  ***
209      !! ** Purpose: Set coupling fields and pass to CICE
210      !!---------------------------------------------------------------------
211      INTEGER, INTENT(in   ) ::   kt   ! ocean time step
212      INTEGER, INTENT(in   ) ::   nsbc ! surface forcing type
213
214      INTEGER  ::   ji, jj, jl                   ! dummy loop indices     
215      REAL(wp), DIMENSION(:,:), POINTER :: ztmp
216      REAL(wp), DIMENSION(:,:,:), POINTER :: ztmpn
217      !!---------------------------------------------------------------------
218
219      IF( nn_timing == 1 )  CALL timing_start('cice_sbc_in')
220      !
221      CALL wrk_alloc( jpi,jpj, ztmp )
222      CALL wrk_alloc( jpi,jpj,ncat, ztmpn )
223
224      IF( kt == nit000 )  THEN
225         IF(lwp) WRITE(numout,*)'cice_sbc_in'
226      ENDIF
227
228      ztmp(:,:)=0.0
229
230! Aggregate ice concentration already set in cice_sbc_out (or cice_sbc_init on
231! the first time-step)
232
233! forced and coupled case
234
235      IF ( (nsbc == 2).OR.(nsbc == 5) ) THEN
236
237         ztmpn(:,:,:)=0.0
238
239! x comp of wind stress (CI_1)
240! U point to F point
241         DO jj=1,jpjm1
242            DO ji=1,jpi
243               ztmp(ji,jj) = 0.5 * (  fr_iu(ji,jj) * utau(ji,jj)      &
244                                    + fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1)
245            ENDDO
246         ENDDO
247         CALL nemo2cice(ztmp,strax,'F', -1. )
248
249! y comp of wind stress (CI_2)
250! V point to F point
251         DO jj=1,jpj
252            DO ji=1,jpim1
253               ztmp(ji,jj) = 0.5 * (  fr_iv(ji,jj) * vtau(ji,jj)      &
254                                    + fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1)
255            ENDDO
256         ENDDO
257         CALL nemo2cice(ztmp,stray,'F', -1. )
258
259! Surface downward latent heat flux (CI_5)
260         IF (nsbc == 2) THEN
261            DO jl=1,ncat
262               ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl)
263            ENDDO
264         ELSE
265! emp_ice is set in sbc_cpl_ice_flx as sublimation-snow
266            qla_ice(:,:,1)= - ( emp_ice(:,:)+sprecip(:,:) ) * Lsub
267! End of temporary code
268            DO jj=1,jpj
269               DO ji=1,jpi
270                  IF (fr_i(ji,jj).eq.0.0) THEN
271                     DO jl=1,ncat
272                        ztmpn(ji,jj,jl)=0.0
273                     ENDDO
274                     ! This will then be conserved in CICE
275                     ztmpn(ji,jj,1)=qla_ice(ji,jj,1)
276                  ELSE
277                     DO jl=1,ncat
278                        ztmpn(ji,jj,jl)=qla_ice(ji,jj,1)*a_i(ji,jj,jl)/fr_i(ji,jj)
279                     ENDDO
280                  ENDIF
281               ENDDO
282            ENDDO
283         ENDIF
284         DO jl=1,ncat
285            CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. )
286
287! GBM conductive flux through ice (CI_6)
288!  Convert to GBM
289            IF (nsbc == 2) THEN
290               ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl)
291            ELSE
292               ztmp(:,:) = botmelt(:,:,jl)
293            ENDIF
294            CALL nemo2cice(ztmp,fcondtopn_f(:,:,jl,:),'T', 1. )
295
296! GBM surface heat flux (CI_7)
297!  Convert to GBM
298            IF (nsbc == 2) THEN
299               ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl) 
300            ELSE
301               ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))
302            ENDIF
303            CALL nemo2cice(ztmp,fsurfn_f(:,:,jl,:),'T', 1. )
304         ENDDO
305
306      ELSE IF (nsbc == 4) THEN
307
308! Pass CORE forcing fields to CICE (which will calculate heat fluxes etc itself)
309! x comp and y comp of atmosphere surface wind (CICE expects on T points)
310         ztmp(:,:) = wndi_ice(:,:)
311         CALL nemo2cice(ztmp,uatm,'T', -1. )
312         ztmp(:,:) = wndj_ice(:,:)
313         CALL nemo2cice(ztmp,vatm,'T', -1. )
314         ztmp(:,:) = SQRT ( wndi_ice(:,:)**2 + wndj_ice(:,:)**2 )
315         CALL nemo2cice(ztmp,wind,'T', 1. )    ! Wind speed (m/s)
316         ztmp(:,:) = qsr_ice(:,:,1)
317         CALL nemo2cice(ztmp,fsw,'T', 1. )     ! Incoming short-wave (W/m^2)
318         ztmp(:,:) = qlw_ice(:,:,1)
319         CALL nemo2cice(ztmp,flw,'T', 1. )     ! Incoming long-wave (W/m^2)
320         ztmp(:,:) = tatm_ice(:,:)
321         CALL nemo2cice(ztmp,Tair,'T', 1. )    ! Air temperature (K)
322         CALL nemo2cice(ztmp,potT,'T', 1. )    ! Potential temp (K)
323! Following line uses MAX(....) to avoid problems if tatm_ice has unset halo rows 
324         ztmp(:,:) = 101000. / ( 287.04 * MAX(1.0,tatm_ice(:,:)) )   
325                                               ! Constant (101000.) atm pressure assumed
326         CALL nemo2cice(ztmp,rhoa,'T', 1. )    ! Air density (kg/m^3)
327         ztmp(:,:) = qatm_ice(:,:)
328         CALL nemo2cice(ztmp,Qa,'T', 1. )      ! Specific humidity (kg/kg)
329         ztmp(:,:)=10.0
330         CALL nemo2cice(ztmp,zlvl,'T', 1. )    ! Atmos level height (m)
331
332! May want to check all values are physically realistic (as in CICE routine
333! prepare_forcing)?
334
335! Divide shortwave into spectral bands (as in prepare_forcing)
336         ztmp(:,:)=qsr_ice(:,:,1)*frcvdr       ! visible direct
337         CALL nemo2cice(ztmp,swvdr,'T', 1. )             
338         ztmp(:,:)=qsr_ice(:,:,1)*frcvdf       ! visible diffuse
339         CALL nemo2cice(ztmp,swvdf,'T', 1. )             
340         ztmp(:,:)=qsr_ice(:,:,1)*frcidr       ! near IR direct
341         CALL nemo2cice(ztmp,swidr,'T', 1. )
342         ztmp(:,:)=qsr_ice(:,:,1)*frcidf       ! near IR diffuse
343         CALL nemo2cice(ztmp,swidf,'T', 1. )
344
345      ENDIF
346
347! Snowfall
348! Ensure fsnow is positive (as in CICE routine prepare_forcing) 
349      ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0) 
350      CALL nemo2cice(ztmp,fsnow,'T', 1. ) 
351
352! Rainfall
353      ztmp(:,:)=fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
354      CALL nemo2cice(ztmp,frain,'T', 1. ) 
355
356! Freezing/melting potential
357! Calculated over NEMO leapfrog timestep (hence 2*dt)
358      nfrzmlt(:,:)=rau0*rcp*fse3t_m(:,:,1)*(Tocnfrz-sst_m(:,:))/(2.0*dt)
359
360      ztmp(:,:) = nfrzmlt(:,:)
361      CALL nemo2cice(ztmp,frzmlt,'T', 1. )
362
363! SST  and SSS
364
365      CALL nemo2cice(sst_m,sst,'T', 1. )
366      CALL nemo2cice(sss_m,sss,'T', 1. )
367
368! x comp and y comp of surface ocean current
369! U point to F point
370      DO jj=1,jpjm1
371         DO ji=1,jpi
372            ztmp(ji,jj)=0.5*(ssu_m(ji,jj)+ssu_m(ji,jj+1))*fmask(ji,jj,1)
373         ENDDO
374      ENDDO
375      CALL nemo2cice(ztmp,uocn,'F', -1. )
376
377! V point to F point
378      DO jj=1,jpj
379         DO ji=1,jpim1
380            ztmp(ji,jj)=0.5*(ssv_m(ji,jj)+ssv_m(ji+1,jj))*fmask(ji,jj,1)
381         ENDDO
382      ENDDO
383      CALL nemo2cice(ztmp,vocn,'F', -1. )
384
385! x comp and y comp of sea surface slope (on F points)
386! T point to F point
387      DO jj=1,jpjm1
388         DO ji=1,jpim1
389            ztmp(ji,jj)=0.5 * (  (ssh_m(ji+1,jj  )-ssh_m(ji,jj  ))/e1u(ji,jj  )   &
390                               + (ssh_m(ji+1,jj+1)-ssh_m(ji,jj+1))/e1u(ji,jj+1) ) & 
391                            *  fmask(ji,jj,1)
392         ENDDO
393      ENDDO
394      CALL nemo2cice(ztmp,ss_tltx,'F', -1. )
395
396! T point to F point
397      DO jj=1,jpjm1
398         DO ji=1,jpim1
399            ztmp(ji,jj)=0.5 * (  (ssh_m(ji  ,jj+1)-ssh_m(ji  ,jj))/e2v(ji  ,jj)   &
400                               + (ssh_m(ji+1,jj+1)-ssh_m(ji+1,jj))/e2v(ji+1,jj) ) &
401                            *  fmask(ji,jj,1)
402         ENDDO
403      ENDDO
404      CALL nemo2cice(ztmp,ss_tlty,'F', -1. )
405
406      CALL wrk_dealloc( jpi,jpj, ztmp )
407      CALL wrk_dealloc( jpi,jpj,ncat, ztmpn )
408      !
409      IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_in')
410      !
411   END SUBROUTINE cice_sbc_in
412
413
414   SUBROUTINE cice_sbc_out (kt,nsbc)
415      !!---------------------------------------------------------------------
416      !!                    ***  ROUTINE cice_sbc_out  ***
417      !! ** Purpose: Get fields from CICE and set surface fields for NEMO
418      !!---------------------------------------------------------------------
419      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
420      INTEGER, INTENT( in  ) ::   nsbc ! surface forcing type
421     
422      INTEGER  ::   ji, jj, jl                 ! dummy loop indices
423      REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
424      !!---------------------------------------------------------------------
425
426      IF( nn_timing == 1 )  CALL timing_start('cice_sbc_out')
427      !
428      CALL wrk_alloc( jpi,jpj, ztmp1, ztmp2 )
429     
430      IF( kt == nit000 )  THEN
431         IF(lwp) WRITE(numout,*)'cice_sbc_out'
432      ENDIF
433     
434! x comp of ocean-ice stress
435      CALL cice2nemo(strocnx,ztmp1,'F', -1. )
436      ss_iou(:,:)=0.0
437! F point to U point
438      DO jj=2,jpjm1
439         DO ji=2,jpim1
440            ss_iou(ji,jj) = 0.5 * ( ztmp1(ji,jj-1) + ztmp1(ji,jj) ) * umask(ji,jj,1)
441         ENDDO
442      ENDDO
443      CALL lbc_lnk( ss_iou , 'U', -1. )
444
445! y comp of ocean-ice stress
446      CALL cice2nemo(strocny,ztmp1,'F', -1. )
447      ss_iov(:,:)=0.0
448! F point to V point
449
450      DO jj=1,jpjm1
451         DO ji=2,jpim1
452            ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1)
453         ENDDO
454      ENDDO
455      CALL lbc_lnk( ss_iov , 'V', -1. )
456
457! x and y comps of surface stress
458! Combine wind stress and ocean-ice stress
459! [Note that fr_iu hasn't yet been updated, so still from start of CICE timestep]
460
461      utau(:,:)=(1.0-fr_iu(:,:))*utau(:,:)-ss_iou(:,:)
462      vtau(:,:)=(1.0-fr_iv(:,:))*vtau(:,:)-ss_iov(:,:)     
463
464! Freshwater fluxes
465
466      IF (nsbc == 2) THEN
467! Note that emp from the forcing files is evap*(1-aice)-(tprecip-aice*sprecip)
468! What we want here is evap*(1-aice)-tprecip*(1-aice) hence manipulation below
469! Not ideal since aice won't be the same as in the atmosphere. 
470! Better to use evap and tprecip? (but for now don't read in evap in this case)
471         emp(:,:)  = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:))
472      ELSE IF (nsbc == 4) THEN
473         emp(:,:)  = (1.0-fr_i(:,:))*emp(:,:)       
474      ELSE IF (nsbc ==5) THEN
475! emp_tot is set in sbc_cpl_ice_flx (call from cice_sbc_in above)
476         emp(:,:) = emp_tot(:,:)+tprecip(:,:)*fr_i(:,:) 
477      ENDIF
478
479      CALL cice2nemo(fresh_gbm,ztmp1,'T', 1. )
480      CALL cice2nemo(fsalt_gbm,ztmp2,'T', 1. )
481! Check to avoid unphysical expression when ice is forming (ztmp1 negative)
482! Otherwise we are effectively allowing ice of higher salinity than the ocean to form
483! which has to be compensated for by the ocean salinity potentially going negative
484! This check breaks conservation but seems reasonable until we have prognostic ice salinity
485! The lines below ensure that when ice is forming emps will not be negative
486! Note the 1000.0 below is to convert from kg salt to g salt (needed for PSU)
487      emps(:,:)=0.0
488      WHERE (ztmp1(:,:).lt.0.0) ztmp2(:,:)=MAX(ztmp2(:,:),ztmp1(:,:)*sss_m(:,:)/1000.0)
489      WHERE (sss_m(:,:).gt.0.0) emps(:,:)=-ztmp1(:,:)+ztmp2(:,:)*1000.0/sss_m(:,:)
490! Now add other non-ice freshwater contributions into emps
491      emps(:,:)=emp(:,:)+emps(:,:)
492 
493      IF (lk_vvl) THEN
494
495! The relevant quantity for calculating the salinity change is emps-emp
496! This is in fact just the ztmp2*1000.0/sss_m term above
497         emp(:,:)=emp(:,:)-ztmp1(:,:)
498
499      ELSE
500
501! Don't remove precip over ice from free surface calculation on basis that the
502! weight of the precip will affect the free surface even if it falls on the ice
503! (same argument that freezing / melting of ice doesn't change the free surface)
504         emp(:,:)  = emp(:,:) - tprecip(:,:)*fr_i(:,:)
505! Sublimation from the ice is treated in a similar way (included in emp but not emps)
506         IF (nsbc == 5 ) THEN
507            emp(:,:) = emp(:,:) + ( emp_ice(:,:) + sprecip(:,:) )
508         ELSE IF (nsbc == 2 ) THEN
509            emp(:,:) = emp(:,:) - qla_ice(:,:,1) / Lsub
510         ENDIF
511
512      ENDIF
513
514      CALL lbc_lnk( emp , 'T', 1. )
515      CALL lbc_lnk( emps , 'T', 1. )
516
517! Solar penetrative radiation and non solar surface heat flux
518
519! Scale qsr and qns according to ice fraction (bulk formulae only)
520
521      IF (nsbc == 4) THEN
522         qsr(:,:)=qsr(:,:)*(1.0-fr_i(:,:))
523         qns(:,:)=qns(:,:)*(1.0-fr_i(:,:))
524      ENDIF
525! Take into account snow melting except for fully coupled when already in qns_tot
526      IF (nsbc == 5) THEN
527         qsr(:,:)= qsr_tot(:,:)
528         qns(:,:)= qns_tot(:,:)
529      ELSE
530         qns(:,:)= qns(:,:)-sprecip(:,:)*Lfresh*(1.0-fr_i(:,:))
531      ENDIF
532
533! Now add in ice / snow related terms
534! [fswthru will be zero unless running with calc_Tsfc=T in CICE]
535      CALL cice2nemo(fswthru_gbm,ztmp1,'T', 1. )
536      qsr(:,:)=qsr(:,:)+ztmp1(:,:)
537      CALL lbc_lnk( qsr , 'T', 1. )
538
539      DO jj=1,jpj
540         DO ji=1,jpi
541            nfrzmlt(ji,jj)=MAX(nfrzmlt(ji,jj),0.0)
542         ENDDO
543      ENDDO
544
545      CALL cice2nemo(fhocn_gbm,ztmp1,'T', 1. )
546      qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:)
547
548      CALL lbc_lnk( qns , 'T', 1. )
549
550! Prepare for the following CICE time-step
551
552      CALL cice2nemo(aice,fr_i,'T', 1. )
553      IF ( (nsbc == 2).OR.(nsbc == 5) ) THEN
554         DO jl=1,ncat
555            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. )
556         ENDDO
557      ENDIF
558
559! T point to U point
560! T point to V point
561      DO jj=1,jpjm1
562         DO ji=1,jpim1
563            fr_iu(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji+1,jj))*umask(ji,jj,1)
564            fr_iv(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji,jj+1))*vmask(ji,jj,1)
565         ENDDO
566      ENDDO
567
568      CALL lbc_lnk ( fr_iu , 'U', 1. )
569      CALL lbc_lnk ( fr_iv , 'V', 1. )
570
571! Release work space
572
573      CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
574      !
575      IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_out')
576      !
577   END SUBROUTINE cice_sbc_out
578
579
580#if defined key_oasis3 || defined key_oasis4
581   SUBROUTINE cice_sbc_hadgam( kt )
582      !!---------------------------------------------------------------------
583      !!                    ***  ROUTINE cice_sbc_hadgam  ***
584      !! ** Purpose: Prepare fields needed to pass to HadGAM3 atmosphere
585      !!
586      !!
587      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
588      !!---------------------------------------------------------------------
589
590      INTEGER  ::   jl                        ! dummy loop index
591      INTEGER  ::   ierror
592
593      IF( nn_timing == 1 )  CALL timing_start('cice_sbc_hadgam')
594      !
595      IF( kt == nit000 )  THEN
596         IF(lwp) WRITE(numout,*)'cice_sbc_hadgam'
597         IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
598      ENDIF
599
600      !                                         ! =========================== !
601      !                                         !   Prepare Coupling fields   !
602      !                                         ! =========================== !
603
604! x and y comp of ice velocity
605
606      CALL cice2nemo(uvel,u_ice,'F', -1. )
607      CALL cice2nemo(vvel,v_ice,'F', -1. )
608
609! Ice concentration (CO_1) = a_i calculated at end of cice_sbc_out 
610
611! Snow and ice thicknesses (CO_2 and CO_3)
612
613      DO jl = 1,ncat
614         CALL cice2nemo(vsnon(:,:,jl,:),ht_s(:,:,jl),'T', 1. )
615         CALL cice2nemo(vicen(:,:,jl,:),ht_i(:,:,jl),'T', 1. )
616      ENDDO
617      !
618      IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_hadgam')
619      !
620   END SUBROUTINE cice_sbc_hadgam
621
622#else
623   SUBROUTINE cice_sbc_hadgam( kt )    ! Dummy routine
624      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
625      WRITE(*,*) 'cice_sbc_hadgam: You should not have seen this print! error?'
626   END SUBROUTINE cice_sbc_hadgam
627#endif
628
629   SUBROUTINE cice_sbc_final
630      !!---------------------------------------------------------------------
631      !!                    ***  ROUTINE cice_sbc_final  ***
632      !! ** Purpose: Finalize CICE
633      !!---------------------------------------------------------------------
634
635      IF(lwp) WRITE(numout,*)'cice_sbc_final'
636
637      CALL CICE_Finalize
638
639   END SUBROUTINE cice_sbc_final
640
641   SUBROUTINE cice_sbc_force (kt)
642      !!---------------------------------------------------------------------
643      !!                    ***  ROUTINE cice_sbc_force  ***
644      !! ** Purpose : Provide CICE forcing from files
645      !!
646      !!---------------------------------------------------------------------
647      !! ** Method  :   READ monthly flux file in NetCDF files
648      !!     
649      !!  snowfall   
650      !!  rainfall   
651      !!  sublimation rate   
652      !!  topmelt (category)
653      !!  botmelt (category)
654      !!
655      !! History :
656      !!----------------------------------------------------------------------
657      !! * Modules used
658      USE iom
659
660      !! * arguments
661      INTEGER, INTENT( in  ) ::   kt ! ocean time step
662
663      INTEGER  ::   ierror             ! return error code
664      INTEGER  ::   ifpr               ! dummy loop index
665      !!
666      CHARACTER(len=100) ::  cn_dir                            !   Root directory for location of CICE forcing files
667      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read
668      TYPE(FLD_N) ::   sn_snow, sn_rain, sn_sblm               ! informations about the fields to be read
669      TYPE(FLD_N) ::   sn_top1, sn_top2, sn_top3, sn_top4, sn_top5
670      TYPE(FLD_N) ::   sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5 
671
672      !!
673      NAMELIST/namsbc_cice/ cn_dir, sn_snow, sn_rain, sn_sblm,   &
674         &                          sn_top1, sn_top2, sn_top3, sn_top4, sn_top5,  &
675         &                          sn_bot1, sn_bot2, sn_bot3, sn_bot4, sn_bot5
676      !!---------------------------------------------------------------------
677
678      !                                         ! ====================== !
679      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
680         !                                      ! ====================== !
681         ! set file information (default values)
682         cn_dir = './'       ! directory in which the model is executed
683
684         ! (NB: frequency positive => hours, negative => months)
685         !            !    file          ! frequency !  variable    ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   !
686         !            !    name          !  (hours)  !   name       !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      !
687         sn_snow = FLD_N( 'snowfall_1m'  ,    -1.    ,  'snowfall'  ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ) 
688         sn_rain = FLD_N( 'rainfall_1m'  ,    -1.    ,  'rainfall'  ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         ) 
689         sn_sblm = FLD_N( 'sublim_1m'    ,    -1.    ,  'sublim'    ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         )
690         sn_top1 = FLD_N( 'topmeltn1_1m' ,    -1.    ,  'topmeltn1' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         )
691         sn_top2 = FLD_N( 'topmeltn2_1m' ,    -1.    ,  'topmeltn2' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         )
692         sn_top3 = FLD_N( 'topmeltn3_1m' ,    -1.    ,  'topmeltn3' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         )
693         sn_top4 = FLD_N( 'topmeltn4_1m' ,    -1.    ,  'topmeltn4' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         )
694         sn_top5 = FLD_N( 'topmeltn5_1m' ,    -1.    ,  'topmeltn5' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         )
695         sn_bot1 = FLD_N( 'botmeltn1_1m' ,    -1.    ,  'botmeltn1' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         )
696         sn_bot2 = FLD_N( 'botmeltn2_1m' ,    -1.    ,  'botmeltn2' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         )
697         sn_bot3 = FLD_N( 'botmeltn3_1m' ,    -1.    ,  'botmeltn3' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         )
698         sn_bot4 = FLD_N( 'botmeltn4_1m' ,    -1.    ,  'botmeltn4' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         )
699         sn_bot5 = FLD_N( 'botmeltn5_1m' ,    -1.    ,  'botmeltn5' ,  .true.    , .true.  ,  ' yearly'  , ''       , ''         )
700
701!         REWIND ( numnam )               ! ... at some point might read in from NEMO namelist?
702!         READ   ( numnam, namsbc_cice )
703
704         ! store namelist information in an array
705         slf_i(jp_snow) = sn_snow   ;   slf_i(jp_rain) = sn_rain   ;   slf_i(jp_sblm) = sn_sblm
706         slf_i(jp_top1) = sn_top1   ;   slf_i(jp_top2) = sn_top2   ;   slf_i(jp_top3) = sn_top3
707         slf_i(jp_top4) = sn_top4   ;   slf_i(jp_top5) = sn_top5   ;   slf_i(jp_bot1) = sn_bot1
708         slf_i(jp_bot2) = sn_bot2   ;   slf_i(jp_bot3) = sn_bot3   ;   slf_i(jp_bot4) = sn_bot4
709         slf_i(jp_bot5) = sn_bot5
710         
711         ! set sf structure
712         ALLOCATE( sf(jpfld), STAT=ierror )
713         IF( ierror > 0 ) THEN
714            CALL ctl_stop( 'cice_sbc_force: unable to allocate sf structure' )   ;   RETURN
715         ENDIF
716
717         DO ifpr= 1, jpfld
718            ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
719            ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
720         END DO
721
722         ! fill sf with slf_i and control print
723         CALL fld_fill( sf, slf_i, cn_dir, 'cice_sbc_force', 'flux formulation for CICE', 'namsbc_cice' )
724         !
725      ENDIF
726
727      CALL fld_read( kt, nn_fsbc, sf )           ! Read input fields and provides the
728      !                                          ! input fields at the current time-step
729
730      ! set the fluxes from read fields
731      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1)
732      tprecip(:,:) = sf(jp_snow)%fnow(:,:,1)+sf(jp_rain)%fnow(:,:,1)
733! May be better to do this conversion somewhere else
734      qla_ice(:,:,1) = -Lsub*sf(jp_sblm)%fnow(:,:,1)
735      topmelt(:,:,1) = sf(jp_top1)%fnow(:,:,1)
736      topmelt(:,:,2) = sf(jp_top2)%fnow(:,:,1)
737      topmelt(:,:,3) = sf(jp_top3)%fnow(:,:,1)
738      topmelt(:,:,4) = sf(jp_top4)%fnow(:,:,1)
739      topmelt(:,:,5) = sf(jp_top5)%fnow(:,:,1)
740      botmelt(:,:,1) = sf(jp_bot1)%fnow(:,:,1)
741      botmelt(:,:,2) = sf(jp_bot2)%fnow(:,:,1)
742      botmelt(:,:,3) = sf(jp_bot3)%fnow(:,:,1)
743      botmelt(:,:,4) = sf(jp_bot4)%fnow(:,:,1)
744      botmelt(:,:,5) = sf(jp_bot5)%fnow(:,:,1)
745
746      ! control print (if less than 100 time-step asked)
747      IF( nitend-nit000 <= 100 .AND. lwp ) THEN
748         WRITE(numout,*) 
749         WRITE(numout,*) '        read forcing fluxes for CICE OK'
750         CALL FLUSH(numout)
751      ENDIF
752
753   END SUBROUTINE cice_sbc_force
754
755   SUBROUTINE nemo2cice( pn, pc, cd_type, psgn)
756      !!---------------------------------------------------------------------
757      !!                    ***  ROUTINE nemo2cice  ***
758      !! ** Purpose :   Transfer field in NEMO array to field in CICE array. 
759#if defined key_nemocice_decomp
760      !!             
761      !!                NEMO and CICE PE sub domains are identical, hence
762      !!                there is no need to gather or scatter data from
763      !!                one PE configuration to another.
764#else
765      !!                Automatically gather/scatter between
766      !!                different processors and blocks
767      !! ** Method :    A. Ensure all haloes are filled in NEMO field (pn)
768      !!                B. Gather pn into global array (png)
769      !!                C. Map png into CICE global array (pcg)
770      !!                D. Scatter pcg to CICE blocks (pc) + update haloes 
771#endif
772      !!---------------------------------------------------------------------
773
774      CHARACTER(len=1), INTENT( in ) ::   &
775          cd_type       ! nature of pn grid-point
776          !             !   = T or F gridpoints
777      REAL(wp), INTENT( in ) ::   &
778          psgn          ! control of the sign change
779          !             !   =-1 , the sign is modified following the type of b.c. used
780          !             !   = 1 , no sign change
781      REAL(wp), DIMENSION(jpi,jpj) :: pn
782#if !defined key_nemocice_decomp
783      REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
784#endif
785      REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
786      INTEGER (int_kind) :: &
787         field_type,        &! id for type of field (scalar, vector, angle)
788         grid_loc            ! id for location on horizontal grid
789                            !  (center, NEcorner, Nface, Eface)
790
791      INTEGER  ::   ji, jj, jn                      ! dummy loop indices
792
793!     A. Ensure all haloes are filled in NEMO field (pn)
794
795      CALL lbc_lnk( pn , cd_type, psgn )
796
797#if defined key_nemocice_decomp
798
799      ! Copy local domain data from NEMO to CICE field
800      pc(:,:,1)=0.0
801      DO jj=2,ny_block
802         DO ji=2,nx_block
803            pc(ji,jj,1)=pn(ji,jj-1)
804         ENDDO
805      ENDDO
806
807#else
808
809!     B. Gather pn into global array (png)
810
811      IF ( jpnij > 1) THEN
812         CALL mppsync
813         CALL mppgather (pn,0,png) 
814         CALL mppsync
815      ELSE
816         png(:,:,1)=pn(:,:)
817      ENDIF
818
819!     C. Map png into CICE global array (pcg)
820
821! Need to make sure this is robust to changes in NEMO halo rows....
822! (may be OK but not 100% sure)
823
824      IF (nproc==0) THEN     
825!        pcg(:,:)=0.0
826         DO jn=1,jpnij
827            DO jj=1,nlcjt(jn)-1
828               DO ji=2,nlcit(jn)-1
829                  pcg(ji+nimppt(jn)-2,jj+njmppt(jn)-1)=png(ji,jj,jn)       
830               ENDDO
831            ENDDO
832         ENDDO
833      ENDIF
834
835#endif
836
837      SELECT CASE ( cd_type )
838         CASE ( 'T' )
839            grid_loc=field_loc_center
840         CASE ( 'F' )                             
841            grid_loc=field_loc_NEcorner
842      END SELECT
843
844      SELECT CASE ( NINT(psgn) )
845         CASE ( -1 )
846            field_type=field_type_vector
847         CASE ( 1 )                             
848            field_type=field_type_scalar
849      END SELECT
850
851#if defined key_nemocice_decomp
852      ! Ensure CICE halos are up to date
853      CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
854#else
855!     D. Scatter pcg to CICE blocks (pc) + update halos
856      CALL scatter_global(pc, pcg, 0, distrb_info, grid_loc, field_type)
857#endif
858
859   END SUBROUTINE nemo2cice
860
861   SUBROUTINE cice2nemo ( pc, pn, cd_type, psgn )
862      !!---------------------------------------------------------------------
863      !!                    ***  ROUTINE cice2nemo  ***
864      !! ** Purpose :   Transfer field in CICE array to field in NEMO array.
865#if defined key_nemocice_decomp
866      !!             
867      !!                NEMO and CICE PE sub domains are identical, hence
868      !!                there is no need to gather or scatter data from
869      !!                one PE configuration to another.
870#else 
871      !!                Automatically deal with scatter/gather between
872      !!                different processors and blocks
873      !! ** Method :    A. Gather CICE blocks (pc) into global array (pcg)
874      !!                B. Map pcg into NEMO global array (png)
875      !!                C. Scatter png into NEMO field (pn) for each processor
876      !!                D. Ensure all haloes are filled in pn
877#endif
878      !!---------------------------------------------------------------------
879
880      CHARACTER(len=1), INTENT( in ) ::   &
881          cd_type       ! nature of pn grid-point
882          !             !   = T or F gridpoints
883      REAL(wp), INTENT( in ) ::   &
884          psgn          ! control of the sign change
885          !             !   =-1 , the sign is modified following the type of b.c. used
886          !             !   = 1 , no sign change
887      REAL(wp), DIMENSION(jpi,jpj) :: pn
888
889#if defined key_nemocice_decomp
890      INTEGER (int_kind) :: &
891         field_type,        & ! id for type of field (scalar, vector, angle)
892         grid_loc             ! id for location on horizontal grid
893                              ! (center, NEcorner, Nface, Eface)
894#else
895      REAL (kind=dbl_kind), dimension(nx_global,ny_global) :: pcg
896#endif
897
898      REAL (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: pc
899
900      INTEGER  ::   ji, jj, jn                      ! dummy loop indices
901
902
903#if defined key_nemocice_decomp
904
905      SELECT CASE ( cd_type )
906         CASE ( 'T' )
907            grid_loc=field_loc_center
908         CASE ( 'F' )                             
909            grid_loc=field_loc_NEcorner
910      END SELECT
911
912      SELECT CASE ( NINT(psgn) )
913         CASE ( -1 )
914            field_type=field_type_vector
915         CASE ( 1 )                             
916            field_type=field_type_scalar
917      END SELECT
918
919      CALL ice_HaloUpdate (pc, halo_info, grid_loc, field_type)
920
921
922      pn(:,:)=0.0
923      DO jj=1,jpjm1
924         DO ji=1,jpim1
925            pn(ji,jj)=pc(ji,jj+1,1)
926         ENDDO
927      ENDDO
928
929#else
930
931!      A. Gather CICE blocks (pc) into global array (pcg)
932
933      CALL gather_global(pcg, pc, 0, distrb_info)
934
935!     B. Map pcg into NEMO global array (png)
936
937! Need to make sure this is robust to changes in NEMO halo rows....
938! (may be OK but not spent much time thinking about it)
939
940      IF (nproc==0) THEN
941         png(:,:,:)=0.0
942         DO jn=1,jpnij
943            DO jj=1,nlcjt(jn)-1
944               DO ji=2,nlcit(jn)-1
945                  png(ji,jj,jn)=pcg(ji+nimppt(jn)-2,jj+njmppt(jn)-1)     
946               ENDDO
947            ENDDO
948         ENDDO
949      ENDIF
950
951!     C. Scatter png into NEMO field (pn) for each processor
952
953      IF ( jpnij > 1) THEN
954         CALL mppsync
955         CALL mppscatter (png,0,pn) 
956         CALL mppsync
957      ELSE
958         pn(:,:)=png(:,:,1)
959      ENDIF
960
961#endif
962
963!     D. Ensure all haloes are filled in pn
964
965      CALL lbc_lnk( pn , cd_type, psgn )
966
967   END SUBROUTINE cice2nemo
968
969#else
970   !!----------------------------------------------------------------------
971   !!   Default option           Dummy module         NO CICE sea-ice model
972   !!----------------------------------------------------------------------
973CONTAINS
974
975   SUBROUTINE sbc_ice_cice ( kt, nsbc )     ! Dummy routine
976      WRITE(*,*) 'sbc_ice_cice: You should not have seen this print! error?', kt
977   END SUBROUTINE sbc_ice_cice
978
979   SUBROUTINE cice_sbc_init (nsbc)    ! Dummy routine
980      WRITE(*,*) 'cice_sbc_init: You should not have seen this print! error?'
981   END SUBROUTINE cice_sbc_init
982
983   SUBROUTINE cice_sbc_final     ! Dummy routine
984      WRITE(*,*) 'cice_sbc_final: You should not have seen this print! error?'
985   END SUBROUTINE cice_sbc_final
986
987#endif
988
989   !!======================================================================
990END MODULE sbcice_cice
Note: See TracBrowser for help on using the repository browser.