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/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 @ 3186

Last change on this file since 3186 was 3186, checked in by smasson, 13 years ago

dev_NEMO_MERGE_2011: replace the old wrk_nemo with the new wrk_nemo

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