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.
Changeset 6363 – NEMO

Changeset 6363


Ignore:
Timestamp:
2016-03-07T16:06:02+01:00 (7 years ago)
Author:
dancopsey
Message:

Merged in dev_r5518_fix_global_ice (which already includes all of dev_r5518_coupling_GSI7_GSI8_landice).

Location:
branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice_bitcomp/NEMOGCM
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice_bitcomp/NEMOGCM/CONFIG/SHARED/field_def.xml

    r5517 r6363  
    193193 
    194194         <!-- * variable related to ice shelf forcing * --> 
     195         <field id="berg_calve"   long_name="Iceberg calving"                               unit="kg/m2/s"  /> 
    195196         <field id="fwfisf"       long_name="Ice shelf melting"                             unit="kg/m2/s"  /> 
    196197         <field id="qisf"         long_name="Ice Shelf Heat Flux"                           unit="W/m2"     /> 
  • branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice_bitcomp/NEMOGCM/CONFIG/SHARED/namelist_ref

    r5501 r6363  
    378378   ln_usecplmask = .false.   !  use a coupling mask file to merge data received from several models 
    379379                             !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel) 
     380   ln_coupled_iceshelf_fluxes = .false. ! If true use rate of change of mass of Greenland and Antarctic icesheets to set the  
     381                                        ! combined magnitude of the iceberg calving and iceshelf melting freshwater fluxes. 
     382   rn_greenland_calving_fraction = 0.5  ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting. 
     383   rn_antarctica_calving_fraction = 0.5 ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting. 
     384   rn_iceshelf_fluxes_tolerance = 1e-6  ! Fractional threshold for detecting differences in icesheet masses (must be positive definite). 
    380385/ 
    381386!----------------------------------------------------------------------- 
  • branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice_bitcomp/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90

    r6362 r6363  
    6868   REAL(wp), PUBLIC ::   rcdsn    =    0.31_wp       !: thermal conductivity of snow 
    6969   REAL(wp), PUBLIC ::   cpic     = 2067.0_wp        !: specific heat for ice  
     70#if defined key_cice 
     71   REAL(wp), PUBLIC ::   lsub     =    2.835e+6_wp   !: pure ice latent heat of sublimation                   [J/kg] 
     72#else 
    7073   REAL(wp), PUBLIC ::   lsub     =    2.834e+6_wp   !: pure ice latent heat of sublimation                   [J/kg] 
     74#endif 
    7175   REAL(wp), PUBLIC ::   lfus     =    0.334e+6_wp   !: latent heat of fusion of fresh ice                    [J/kg] 
    7276   REAL(wp), PUBLIC ::   tmut     =    0.054_wp      !: decrease of seawater meltpoint with salinity 
  • branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice_bitcomp/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90

    r6362 r6363  
    2525   USE icbutl         ! iceberg utility routines 
    2626 
     27   USE sbc_oce        ! for icesheet freshwater input variables 
     28   USE in_out_manager 
     29   USE iom 
     30 
    2731   IMPLICIT NONE 
    2832   PRIVATE 
     
    4852      ! 
    4953      REAL(wp)                        :: zcalving_used, zdist, zfact 
     54      REAL(wp)                        :: zgreenland_calving_sum, zantarctica_calving_sum 
    5055      INTEGER                         :: jn, ji, jj                    ! loop counters 
    5156      INTEGER                         :: imx                           ! temporary integer for max berg class 
     
    5964      zfact = ( (1000._wp)**3 / ( NINT(rday) * nyear_len(1) ) ) * 850._wp 
    6065      berg_grid%calving(:,:) = src_calving(:,:) * tmask_i(:,:) * zfact 
     66 
     67      IF( lk_oasis) THEN 
     68      ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true 
     69      IF( ln_coupled_iceshelf_fluxes ) THEN 
     70 
     71        ! Adjust total calving rates so that sum of iceberg calving and iceshelf melting in the northern 
     72        ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets 
     73        ! to preserve total freshwater conservation in coupled models without an active ice sheet model. 
     74 
     75         zgreenland_calving_sum = SUM( berg_grid%calving(:,:) * greenland_icesheet_mask(:,:) ) 
     76         IF( lk_mpp ) CALL mpp_sum( zgreenland_calving_sum ) 
     77         WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                                 & 
     78        &    berg_grid%calving(:,:) = berg_grid%calving(:,:) * greenland_icesheet_mass_rate_of_change * rn_greenland_calving_fraction & 
     79        &                                     / ( zgreenland_calving_sum + 1.0e-10_wp ) 
     80 
     81         ! check 
     82         IF(lwp) WRITE(numout, *) 'Greenland iceberg calving climatology (kg/s) : ',zgreenland_calving_sum 
     83         zgreenland_calving_sum = SUM( berg_grid%calving(:,:) * greenland_icesheet_mask(:,:) ) 
     84         IF( lk_mpp ) CALL mpp_sum( zgreenland_calving_sum ) 
     85         IF(lwp) WRITE(numout, *) 'Greenland iceberg calving adjusted value (kg/s) : ',zgreenland_calving_sum 
     86 
     87         zantarctica_calving_sum = SUM( berg_grid%calving(:,:) * antarctica_icesheet_mask(:,:) ) 
     88         IF( lk_mpp ) CALL mpp_sum( zantarctica_calving_sum ) 
     89         WHERE( antarctica_icesheet_mask(:,:) == 1.0 )                                                                              & 
     90         berg_grid%calving(:,:) = berg_grid%calving(:,:) * antarctica_icesheet_mass_rate_of_change * rn_antarctica_calving_fraction & 
     91        &                           / ( zantarctica_calving_sum + 1.0e-10_wp ) 
     92  
     93         ! check 
     94         IF(lwp) WRITE(numout, *) 'Antarctica iceberg calving climatology (kg/s) : ',zantarctica_calving_sum 
     95         zantarctica_calving_sum = SUM( berg_grid%calving(:,:) * antarctica_icesheet_mask(:,:) ) 
     96         IF( lk_mpp ) CALL mpp_sum( zantarctica_calving_sum ) 
     97         IF(lwp) WRITE(numout, *) 'Antarctica iceberg calving adjusted value (kg/s) : ',zantarctica_calving_sum 
     98 
     99      ENDIF 
     100      ENDIF 
     101    
     102      CALL iom_put( 'berg_calve', berg_grid%calving(:,:) ) 
    61103 
    62104      ! Heat in units of W/m2, and mask (just in case) 
  • branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice_bitcomp/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r6362 r6363  
    2424   USE trdmxl_oce      ! ocean active mixed layer tracers trends variables 
    2525   USE divcur          ! hor. divergence and curl      (div & cur routines) 
     26   USE sbc_oce         ! for icesheet freshwater input variables 
    2627 
    2728   IMPLICIT NONE 
     
    145146                     CALL iom_rstput( kt, nitrst, numrow, 'rhd'    , rhd       ) 
    146147#endif 
     148                     IF( lk_oasis) THEN 
     149                     ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true 
     150                     IF( ln_coupled_iceshelf_fluxes ) THEN 
     151                        CALL iom_rstput( kt, nitrst, numrow, 'greenland_icesheet_mass', greenland_icesheet_mass ) 
     152                        CALL iom_rstput( kt, nitrst, numrow, 'greenland_icesheet_timelapsed', greenland_icesheet_timelapsed ) 
     153                        CALL iom_rstput( kt, nitrst, numrow, 'greenland_icesheet_mass_roc', greenland_icesheet_mass_rate_of_change ) 
     154                        CALL iom_rstput( kt, nitrst, numrow, 'antarctica_icesheet_mass', antarctica_icesheet_mass ) 
     155                        CALL iom_rstput( kt, nitrst, numrow, 'antarctica_icesheet_timelapsed', antarctica_icesheet_timelapsed ) 
     156                        CALL iom_rstput( kt, nitrst, numrow, 'antarctica_icesheet_mass_roc', antarctica_icesheet_mass_rate_of_change ) 
     157                     ENDIF 
     158                     ENDIF 
     159 
    147160      IF( kt == nitrst ) THEN 
    148161         CALL iom_close( numrow )     ! close the restart file (only at last time step) 
     
    258271#endif 
    259272      ! 
     273      IF( iom_varid( numror, 'greenland_icesheet_mass', ldstop = .FALSE. ) > 0 )   THEN 
     274         CALL iom_get( numror, 'greenland_icesheet_mass', greenland_icesheet_mass ) 
     275         CALL iom_get( numror, 'greenland_icesheet_timelapsed', greenland_icesheet_timelapsed ) 
     276         CALL iom_get( numror, 'greenland_icesheet_mass_roc', greenland_icesheet_mass_rate_of_change ) 
     277      ELSE 
     278         greenland_icesheet_mass = 0.0  
     279         greenland_icesheet_mass_rate_of_change = 0.0  
     280         greenland_icesheet_timelapsed = 0.0 
     281      ENDIF 
     282      IF( iom_varid( numror, 'antarctica_icesheet_mass', ldstop = .FALSE. ) > 0 )   THEN 
     283         CALL iom_get( numror, 'antarctica_icesheet_mass', antarctica_icesheet_mass ) 
     284         CALL iom_get( numror, 'antarctica_icesheet_timelapsed', antarctica_icesheet_timelapsed ) 
     285         CALL iom_get( numror, 'antarctica_icesheet_mass_roc', antarctica_icesheet_mass_rate_of_change ) 
     286      ELSE 
     287         antarctica_icesheet_mass = 0.0  
     288         antarctica_icesheet_mass_rate_of_change = 0.0  
     289         antarctica_icesheet_timelapsed = 0.0 
     290      ENDIF 
     291 
    260292      IF( neuler == 0 ) THEN                                  ! Euler restart (neuler=0) 
    261293         tsb  (:,:,:,:) = tsn  (:,:,:,:)                             ! all before fields set to now values 
  • branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice_bitcomp/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

    r6362 r6363  
    101101   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr_iu              !: ice fraction at NEMO U point 
    102102   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr_iv              !: ice fraction at NEMO V point 
    103     
     103   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sstfrz             !: sea surface freezing temperature 
     104   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tsfc_ice           !: sea-ice surface skin temperature (on categories) 
     105   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   kn_ice             !: sea-ice surface layer thermal conductivity (on cats) 
     106 
    104107   ! variables used in the coupled interface 
    105108   INTEGER , PUBLIC, PARAMETER ::   jpl = ncat 
    106109   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice          ! jpi, jpj 
     110   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  a_p, ht_p ! Meltpond fraction and depth 
    107111#endif 
    108112    
     
    152156 
    153157#if defined key_cice 
    154       ALLOCATE( qla_ice(jpi,jpj,1)    , qlw_ice(jpi,jpj,1)    , qsr_ice(jpi,jpj,1)    , & 
     158      ALLOCATE( qla_ice(jpi,jpj,ncat) , qlw_ice(jpi,jpj,1)    , qsr_ice(jpi,jpj,1)    , & 
    155159                wndi_ice(jpi,jpj)     , tatm_ice(jpi,jpj)     , qatm_ice(jpi,jpj)     , & 
    156160                wndj_ice(jpi,jpj)     , nfrzmlt(jpi,jpj)      , ss_iou(jpi,jpj)       , & 
    157161                ss_iov(jpi,jpj)       , fr_iu(jpi,jpj)        , fr_iv(jpi,jpj)        , & 
    158162                a_i(jpi,jpj,ncat)     , topmelt(jpi,jpj,ncat) , botmelt(jpi,jpj,ncat) , & 
    159                 STAT= ierr(1) ) 
    160       IF( ln_cpl )   ALLOCATE( u_ice(jpi,jpj)        , fr1_i0(jpi,jpj)       , tn_ice (jpi,jpj,1)    , & 
     163                sstfrz(jpi,jpj)       , STAT= ierr(1) ) 
     164   ! Alex West: Allocating tn_ice with 5 categories.  When NEMO is used with CICE, this variable 
     165   ! represents top layer ice temperature, which is multi-category. 
     166      IF( ln_cpl )   ALLOCATE( u_ice(jpi,jpj)        , fr1_i0(jpi,jpj)       , tn_ice (jpi,jpj,jpl)  , & 
    161167         &                     v_ice(jpi,jpj)        , fr2_i0(jpi,jpj)       , alb_ice(jpi,jpj,1)    , & 
    162168         &                     emp_ice(jpi,jpj)      , qns_ice(jpi,jpj,1)    , dqns_ice(jpi,jpj,1)   , & 
    163          &                     STAT= ierr(2) ) 
     169         &                     a_p(jpi,jpj,jpl)      , ht_p(jpi,jpj,jpl)     , tsfc_ice(jpi,jpj,jpl) , & 
     170         &                     kn_ice(jpi,jpj,jpl) ,    STAT=ierr(2) ) 
    164171       
    165172#endif 
  • branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice_bitcomp/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r6362 r6363  
    125125#endif 
    126126   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl) 
     127   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   greenland_icesheet_mass_array, greenland_icesheet_mask 
     128   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   antarctica_icesheet_mass_array, antarctica_icesheet_mask 
    127129 
    128130   !!---------------------------------------------------------------------- 
     
    137139   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3t_m     !: mean (nn_fsbc time-step) sea surface layer thickness       [m] 
    138140   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frq_m     !: mean (nn_fsbc time-step) fraction of solar net radiation absorbed in the 1st T level [-] 
     141    
     142   !!---------------------------------------------------------------------- 
     143   !!  Surface scalars of total ice sheet mass for Greenland and Antarctica,  
     144   !! passed from atmosphere to be converted to dvol and hence a freshwater  
     145   !! flux  by using old values. New values are saved in the dump, to become 
     146   !! old values next coupling timestep. Freshwater fluxes split between  
     147   !! sub iceshelf melting and iceberg calving, scalled to flux per second 
     148   !!---------------------------------------------------------------------- 
     149    
     150   REAL(wp), PUBLIC  :: greenland_icesheet_mass, greenland_icesheet_mass_rate_of_change, greenland_icesheet_timelapsed  
     151   REAL(wp), PUBLIC  :: antarctica_icesheet_mass, antarctica_icesheet_mass_rate_of_change, antarctica_icesheet_timelapsed 
     152 
     153   ! sbccpl namelist parameters associated with icesheet freshwater input code. Included here rather than in sbccpl.F90 to  
     154   ! avoid circular dependencies. 
     155   LOGICAL, PUBLIC     ::   ln_coupled_iceshelf_fluxes     ! If true use rate of change of mass of Greenland and Antarctic icesheets to set the  
     156                                                           ! combined magnitude of the iceberg calving and iceshelf melting freshwater fluxes. 
     157   REAL(wp), PUBLIC    ::   rn_greenland_calving_fraction  ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting. 
     158   REAL(wp), PUBLIC    ::   rn_antarctica_calving_fraction ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting. 
     159   REAL(wp), PUBLIC    ::   rn_iceshelf_fluxes_tolerance   ! Absolute tolerance for detecting differences in icesheet masses.  
    139160 
    140161   !! * Substitutions 
     
    172193         &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) ,      & 
    173194         &      ssv_m  (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) ) 
     195      ALLOCATE( greenland_icesheet_mass_array(jpi,jpj) , antarctica_icesheet_mass_array(jpi,jpj) ) 
     196      ALLOCATE( greenland_icesheet_mask(jpi,jpj) , antarctica_icesheet_mask(jpi,jpj) ) 
    174197         ! 
    175198#if defined key_vvl 
  • branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice_bitcomp/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r6362 r6363  
    4646   USE p4zflx, ONLY : oce_co2 
    4747#endif 
    48 #if defined key_cice 
    49    USE ice_domain_size, only: ncat 
    50 #endif 
    5148#if defined key_lim3 
    5249   USE limthd_dh       ! for CALL lim_thd_snwblow 
    5350#endif 
     51   USE lib_fortran, ONLY: glob_sum 
    5452 
    5553   IMPLICIT NONE 
     
    105103   INTEGER, PARAMETER ::   jpr_e3t1st = 41            ! first T level thickness  
    106104   INTEGER, PARAMETER ::   jpr_fraqsr = 42            ! fraction of solar net radiation absorbed in the first ocean level 
    107    INTEGER, PARAMETER ::   jprcv      = 42            ! total number of fields received 
     105   INTEGER, PARAMETER ::   jpr_ts_ice = 43            ! skin temperature of sea-ice (used for melt-ponds) 
     106   INTEGER, PARAMETER ::   jpr_grnm   = 44            ! Greenland ice mass 
     107   INTEGER, PARAMETER ::   jpr_antm   = 45            ! Antarctic ice mass 
     108   INTEGER, PARAMETER ::   jprcv      = 45            ! total number of fields received 
    108109 
    109110   INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere 
     
    135136   INTEGER, PARAMETER ::   jps_e3t1st = 27            ! first level depth (vvl) 
    136137   INTEGER, PARAMETER ::   jps_fraqsr = 28            ! fraction of solar net radiation absorbed in the first ocean level 
    137    INTEGER, PARAMETER ::   jpsnd      = 28            ! total number of fields sended 
     138   INTEGER, PARAMETER ::   jps_a_p    = 29            ! meltpond fraction   
     139   INTEGER, PARAMETER ::   jps_ht_p   = 30            ! meltpond depth (m)  
     140   INTEGER, PARAMETER ::   jps_kice   = 31            ! ice surface layer thermal conductivity 
     141   INTEGER, PARAMETER ::   jps_sstfrz = 32            ! sea-surface freezing temperature 
     142   INTEGER, PARAMETER ::   jps_fice1  = 33            ! first-order ice concentration (for time-travelling ice coupling) 
     143   INTEGER, PARAMETER ::   jpsnd      = 33            ! total number of fields sended 
    138144 
    139145   !                                                         !!** namelist namsbc_cpl ** 
     
    146152   END TYPE FLD_C 
    147153   ! Send to the atmosphere                           ! 
    148    TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2                         
     154   TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2, sn_snd_cond, sn_snd_mpnd, sn_snd_sstfrz, sn_snd_thick1 
     155 
    149156   ! Received from the atmosphere                     ! 
    150157   TYPE(FLD_C) ::   sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, sn_rcv_qns, sn_rcv_emp, sn_rcv_rnf 
    151    TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2                         
     158   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_ts_ice, sn_rcv_grnm, sn_rcv_antm 
    152159   ! Other namelist parameters                        ! 
    153160   INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data 
     
    216223      REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos 
    217224      !! 
    218       NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2,      & 
    219          &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr,      & 
    220          &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal   , sn_rcv_iceflx,   & 
    221          &                  sn_rcv_co2 , nn_cplmodel  , ln_usecplmask 
     225      NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick , sn_snd_crt   , sn_snd_co2,     & 
     226         &                  sn_snd_cond, sn_snd_mpnd  , sn_snd_sstfrz, sn_snd_thick1,                 & 
     227         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,     & 
     228         &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   , sn_rcv_iceflx,  & 
     229         &                  sn_rcv_co2 , sn_rcv_grnm  , sn_rcv_antm  , sn_rcv_ts_ice, nn_cplmodel  ,  & 
     230         &                  ln_usecplmask, ln_coupled_iceshelf_fluxes, rn_greenland_calving_fraction, & 
     231         &                  rn_antarctica_calving_fraction, rn_iceshelf_fluxes_tolerance 
    222232      !!--------------------------------------------------------------------- 
    223233      ! 
     
    258268         WRITE(numout,*)'      runoffs                         = ', TRIM(sn_rcv_rnf%cldes   ), ' (', TRIM(sn_rcv_rnf%clcat   ), ')' 
    259269         WRITE(numout,*)'      calving                         = ', TRIM(sn_rcv_cal%cldes   ), ' (', TRIM(sn_rcv_cal%clcat   ), ')' 
     270         WRITE(numout,*)'      Greenland ice mass              = ', TRIM(sn_rcv_grnm%cldes  ), ' (', TRIM(sn_rcv_grnm%clcat  ), ')' 
     271         WRITE(numout,*)'      Antarctica ice mass             = ', TRIM(sn_rcv_antm%cldes  ), ' (', TRIM(sn_rcv_antm%clcat  ), ')' 
    260272         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')' 
    261273         WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')' 
     
    269281         WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd 
    270282         WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')' 
     283         WRITE(numout,*)'      ice effective conductivity      = ', TRIM(sn_snd_cond%cldes   ), ' (', TRIM(sn_snd_cond%clcat   ), ')' 
     284         WRITE(numout,*)'      meltponds fraction & depth      = ', TRIM(sn_snd_mpnd%cldes  ), ' (', TRIM(sn_snd_mpnd%clcat   ), ')' 
     285         WRITE(numout,*)'      sea surface freezing temp       = ', TRIM(sn_snd_sstfrz%cldes   ), ' (', TRIM(sn_snd_sstfrz%clcat   ), ')' 
     286 
    271287         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel 
    272288         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask 
     289         WRITE(numout,*)'  ln_coupled_iceshelf_fluxes          = ', ln_coupled_iceshelf_fluxes 
     290         WRITE(numout,*)'  rn_greenland_calving_fraction       = ', rn_greenland_calving_fraction 
     291         WRITE(numout,*)'  rn_antarctica_calving_fraction      = ', rn_antarctica_calving_fraction 
     292         WRITE(numout,*)'  rn_iceshelf_fluxes_tolerance        = ', rn_iceshelf_fluxes_tolerance 
    273293      ENDIF 
    274294 
     
    383403      srcv(jpr_snow)%clname = 'OTotSnow'      ! Snow = solid precipitation 
    384404      srcv(jpr_tevp)%clname = 'OTotEvap'      ! total evaporation (over oce + ice sublimation) 
    385       srcv(jpr_ievp)%clname = 'OIceEvap'      ! evaporation over ice = sublimation 
     405      srcv(jpr_ievp)%clname = 'OIceEvp'      ! evaporation over ice = sublimation 
    386406      srcv(jpr_sbpr)%clname = 'OSubMPre'      ! sublimation - liquid precipitation - solid precipitation  
    387407      srcv(jpr_semp)%clname = 'OISubMSn'      ! ice solid water budget = sublimation - solid precipitation 
     
    396416      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' ) 
    397417      END SELECT 
    398  
     418      !Set the number of categories for coupling of sublimation 
     419      IF ( TRIM( sn_rcv_emp%clcat ) == 'yes' ) srcv(jpr_ievp)%nct = jpl 
     420      ! 
    399421      !                                                      ! ------------------------- ! 
    400422      !                                                      !     Runoffs & Calving     !    
     
    410432      ! 
    411433      srcv(jpr_cal   )%clname = 'OCalving'   ;   IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE. 
     434      srcv(jpr_grnm  )%clname = 'OGrnmass'   ;   IF( TRIM( sn_rcv_grnm%cldes ) == 'coupled' )   srcv(jpr_grnm)%laction = .TRUE. 
     435      srcv(jpr_antm  )%clname = 'OAntmass'   ;   IF( TRIM( sn_rcv_antm%cldes ) == 'coupled' )   srcv(jpr_antm)%laction = .TRUE. 
     436 
    412437 
    413438      !                                                      ! ------------------------- ! 
     
    483508         srcv(jpr_topm:jpr_botm)%laction = .TRUE. 
    484509      ENDIF 
     510       
     511#if defined key_cice && ! defined key_cice4 
     512      !                                                      ! ----------------------------- ! 
     513      !                                                      !  sea-ice skin temperature     !    
     514      !                                                      !  used in meltpond scheme      ! 
     515      !                                                      !  May be calculated in Atm     ! 
     516      !                                                      ! ----------------------------- ! 
     517      srcv(jpr_ts_ice)%clname = 'OTsfIce' 
     518      IF ( TRIM( sn_rcv_ts_ice%cldes ) == 'ice' ) srcv(jpr_ts_ice)%laction = .TRUE. 
     519      IF ( TRIM( sn_rcv_ts_ice%clcat ) == 'yes' ) srcv(jpr_ts_ice)%nct = jpl 
     520      !TODO: Should there be a consistency check here? 
     521#endif 
     522 
    485523      !                                                      ! ------------------------------- ! 
    486524      !                                                      !   OPA-SAS coupling - rcv by opa !    
     
    600638      !                                                      ! ------------------------- ! 
    601639      ssnd(jps_toce)%clname = 'O_SSTSST' 
    602       ssnd(jps_tice)%clname = 'O_TepIce' 
     640      ssnd(jps_tice)%clname = 'OTepIce' 
    603641      ssnd(jps_tmix)%clname = 'O_TepMix' 
    604642      SELECT CASE( TRIM( sn_snd_temp%cldes ) ) 
    605643      CASE( 'none'                                 )       ! nothing to do 
    606644      CASE( 'oce only'                             )   ;   ssnd( jps_toce )%laction = .TRUE. 
    607       CASE( 'oce and ice' , 'weighted oce and ice' ) 
     645      CASE( 'oce and ice' , 'weighted oce and ice' , 'oce and weighted ice') 
    608646         ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE. 
    609647         IF ( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = jpl 
     
    634672 
    635673      !                                                      ! ------------------------- ! 
    636       !                                                      !  Ice fraction & Thickness !  
     674      !                                                      !  Ice fraction & Thickness  
    637675      !                                                      ! ------------------------- ! 
    638676      ssnd(jps_fice)%clname = 'OIceFrc' 
    639677      ssnd(jps_hice)%clname = 'OIceTck' 
    640678      ssnd(jps_hsnw)%clname = 'OSnwTck' 
     679      ssnd(jps_a_p)%clname  = 'OPndFrc' 
     680      ssnd(jps_ht_p)%clname = 'OPndTck' 
     681      ssnd(jps_fice1)%clname = 'OIceFrd' 
    641682      IF( k_ice /= 0 ) THEN 
    642683         ssnd(jps_fice)%laction = .TRUE.                  ! if ice treated in the ocean (even in climato case) 
     684         ssnd(jps_fice1)%laction = .TRUE.                 ! First-order regridded ice concentration, to be used 
     685                                                     ! in producing atmos-to-ice fluxes 
    643686! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now 
    644687         IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_fice)%nct = jpl 
     688         IF ( TRIM( sn_snd_thick1%clcat ) == 'yes' ) ssnd(jps_fice1)%nct = jpl 
    645689      ENDIF 
    646690       
     
    657701      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_thick%cldes' ) 
    658702      END SELECT 
     703 
     704      !                                                      ! ------------------------- ! 
     705      !                                                      ! Ice Meltponds             ! 
     706      !                                                      ! ------------------------- ! 
     707#if defined key_cice && ! defined key_cice4 
     708      ! Meltponds only CICE5  
     709      ssnd(jps_a_p)%clname = 'OPndFrc'    
     710      ssnd(jps_ht_p)%clname = 'OPndTck'    
     711      SELECT CASE ( TRIM( sn_snd_mpnd%cldes ) ) 
     712      CASE ( 'none' ) 
     713         ssnd(jps_a_p)%laction = .FALSE. 
     714         ssnd(jps_ht_p)%laction = .FALSE. 
     715      CASE ( 'ice only' )  
     716         ssnd(jps_a_p)%laction = .TRUE. 
     717         ssnd(jps_ht_p)%laction = .TRUE. 
     718         IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN 
     719            ssnd(jps_a_p)%nct = jpl 
     720            ssnd(jps_ht_p)%nct = jpl 
     721         ELSE 
     722            IF ( jpl > 1 ) THEN 
     723               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_mpnd%cldes if not exchanging category fields' ) 
     724            ENDIF 
     725         ENDIF 
     726      CASE ( 'weighted ice' )  
     727         ssnd(jps_a_p)%laction = .TRUE. 
     728         ssnd(jps_ht_p)%laction = .TRUE. 
     729         IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN 
     730            ssnd(jps_a_p)%nct = jpl  
     731            ssnd(jps_ht_p)%nct = jpl  
     732         ENDIF 
     733      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_mpnd%cldes' ) 
     734      END SELECT 
     735#else 
     736      IF( TRIM( sn_snd_mpnd%cldes /= 'none' ) THEN 
     737         CALL ctl_stop('Meltponds can only be used with CICEv5') 
     738      ENDIF 
     739#endif 
    659740 
    660741      !                                                      ! ------------------------- ! 
     
    689770      !                                                      ! ------------------------- ! 
    690771      ssnd(jps_co2)%clname = 'O_CO2FLX' ;  IF( TRIM(sn_snd_co2%cldes) == 'coupled' )    ssnd(jps_co2 )%laction = .TRUE. 
     772      ! 
     773       
     774      !                                                      ! ------------------------- ! 
     775      !                                                      ! Sea surface freezing temp ! 
     776      !                                                      ! ------------------------- ! 
     777      ssnd(jps_sstfrz)%clname = 'O_SSTFrz' ; IF( TRIM(sn_snd_sstfrz%cldes) == 'coupled' )  ssnd(jps_sstfrz)%laction = .TRUE. 
     778      ! 
     779      !                                                      ! ------------------------- ! 
     780      !                                                      !    Ice conductivity       ! 
     781      !                                                      ! ------------------------- ! 
     782      ! Note that ultimately we will move to passing an ocean effective conductivity as well so there 
     783      ! will be some changes to the parts of the code which currently relate only to ice conductivity 
     784      ssnd(jps_kice )%clname = 'OIceKn' 
     785      SELECT CASE ( TRIM( sn_snd_cond%cldes ) ) 
     786      CASE ( 'none' ) 
     787         ssnd(jps_kice)%laction = .FALSE. 
     788      CASE ( 'ice only' ) 
     789         ssnd(jps_kice)%laction = .TRUE. 
     790         IF ( TRIM( sn_snd_cond%clcat ) == 'yes' ) THEN 
     791            ssnd(jps_kice)%nct = jpl 
     792         ELSE 
     793            IF ( jpl > 1 ) THEN 
     794               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_cond%cldes if not exchanging category fields' ) 
     795            ENDIF 
     796         ENDIF 
     797      CASE ( 'weighted ice' ) 
     798         ssnd(jps_kice)%laction = .TRUE. 
     799         IF ( TRIM( sn_snd_cond%clcat ) == 'yes' ) ssnd(jps_kice)%nct = jpl 
     800      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_cond%cldes' ) 
     801      END SELECT 
     802      ! 
     803       
    691804 
    692805      !                                                      ! ------------------------------- ! 
     
    785898      ncpl_qsr_freq = 86400 / ncpl_qsr_freq 
    786899 
     900      IF( ln_coupled_iceshelf_fluxes ) THEN 
     901          ! Crude masks to separate the Antarctic and Greenland icesheets. Obviously something 
     902          ! more complicated could be done if required. 
     903          greenland_icesheet_mask = 0.0 
     904          WHERE( gphit >= 0.0 ) greenland_icesheet_mask = 1.0 
     905          antarctica_icesheet_mask = 0.0 
     906          WHERE( gphit < 0.0 ) antarctica_icesheet_mask = 1.0 
     907 
     908          ! initialise other variables 
     909          greenland_icesheet_mass_array(:,:) = 0.0 
     910          antarctica_icesheet_mass_array(:,:) = 0.0 
     911 
     912          IF( .not. ln_rstart ) THEN 
     913             greenland_icesheet_mass = 0.0  
     914             greenland_icesheet_mass_rate_of_change = 0.0  
     915             greenland_icesheet_timelapsed = 0.0 
     916             antarctica_icesheet_mass = 0.0  
     917             antarctica_icesheet_mass_rate_of_change = 0.0  
     918             antarctica_icesheet_timelapsed = 0.0 
     919          ENDIF 
     920 
     921      ENDIF 
     922 
    787923      CALL wrk_dealloc( jpi,jpj, zacs, zaos ) 
    788924      ! 
     
    843979      !! 
    844980      LOGICAL  ::   llnewtx, llnewtau      ! update wind stress components and module?? 
    845       INTEGER  ::   ji, jj, jn             ! dummy loop indices 
     981      INTEGER  ::   ji, jj, jl, jn         ! dummy loop indices 
    846982      INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdttra did not change since nit000) 
    847983      REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars      
     984      REAL(wp) ::   zgreenland_icesheet_mass_in, zantarctica_icesheet_mass_in 
     985      REAL(wp) ::   zgreenland_icesheet_mass_b, zantarctica_icesheet_mass_b 
     986      REAL(wp) ::   zmask_sum, zepsilon       
    848987      REAL(wp) ::   zcoef                  ! temporary scalar 
    849988      REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3 
     
    9971136#endif 
    9981137 
     1138#if defined key_cice && ! defined key_cice4 
     1139      !  ! Sea ice surface skin temp: 
     1140      IF( srcv(jpr_ts_ice)%laction ) THEN 
     1141        DO jl = 1, jpl 
     1142          DO jj = 1, jpj 
     1143            DO ji = 1, jpi 
     1144              IF (frcv(jpr_ts_ice)%z3(ji,jj,jl) > 0.0) THEN 
     1145                tsfc_ice(ji,jj,jl) = 0.0 
     1146              ELSE IF (frcv(jpr_ts_ice)%z3(ji,jj,jl) < -60.0) THEN 
     1147                tsfc_ice(ji,jj,jl) = -60.0 
     1148              ELSE 
     1149                tsfc_ice(ji,jj,jl) = frcv(jpr_ts_ice)%z3(ji,jj,jl) 
     1150              ENDIF 
     1151            END DO 
     1152          END DO 
     1153        END DO 
     1154      ENDIF 
     1155#endif 
     1156 
    9991157      !  Fields received by SAS when OASIS coupling 
    10001158      !  (arrays no more filled at sbcssm stage) 
     
    11101268 
    11111269      ENDIF 
     1270       
     1271      !                                                        ! land ice masses : Greenland 
     1272      zepsilon = rn_iceshelf_fluxes_tolerance 
     1273 
     1274 
     1275      ! See if we need zmask_sum... 
     1276      IF ( srcv(jpr_grnm)%laction .OR. srcv(jpr_antm)%laction ) THEN 
     1277         zmask_sum = glob_sum( tmask(:,:,1) ) 
     1278      ENDIF 
     1279 
     1280      IF( srcv(jpr_grnm)%laction ) THEN 
     1281         greenland_icesheet_mass_array(:,:) = frcv(jpr_grnm)%z3(:,:,1) 
     1282         ! take average over ocean points of input array to avoid cumulative error over time 
     1283 
     1284         ! The following must be bit reproducible over different PE decompositions 
     1285         zgreenland_icesheet_mass_in = glob_sum( greenland_icesheet_mass_array(:,:) * tmask(:,:,1) ) 
     1286 
     1287         zgreenland_icesheet_mass_in = zgreenland_icesheet_mass_in / zmask_sum 
     1288         greenland_icesheet_timelapsed = greenland_icesheet_timelapsed + rdt          
     1289         IF( ABS( zgreenland_icesheet_mass_in - greenland_icesheet_mass ) > zepsilon ) THEN 
     1290            zgreenland_icesheet_mass_b = greenland_icesheet_mass 
     1291             
     1292            ! Only update the mass if it has increased 
     1293            IF ( (zgreenland_icesheet_mass_in - greenland_icesheet_mass) > 0.0 ) THEN 
     1294               greenland_icesheet_mass = zgreenland_icesheet_mass_in 
     1295            ENDIF 
     1296             
     1297            IF( zgreenland_icesheet_mass_b /= 0.0 ) & 
     1298           &     greenland_icesheet_mass_rate_of_change = ( greenland_icesheet_mass - zgreenland_icesheet_mass_b ) / greenland_icesheet_timelapsed  
     1299            greenland_icesheet_timelapsed = 0.0_wp        
     1300         ENDIF 
     1301         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass (kg) read in is ', zgreenland_icesheet_mass_in 
     1302         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass (kg) used is    ', greenland_icesheet_mass 
     1303         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass rate of change (kg/s) is ', greenland_icesheet_mass_rate_of_change 
     1304         IF(lwp) WRITE(numout,*) 'Greenland icesheet seconds lapsed since last change is ', greenland_icesheet_timelapsed 
     1305      ENDIF 
     1306 
     1307      !                                                        ! land ice masses : Antarctica 
     1308      IF( srcv(jpr_antm)%laction ) THEN 
     1309         antarctica_icesheet_mass_array(:,:) = frcv(jpr_antm)%z3(:,:,1) 
     1310         ! take average over ocean points of input array to avoid cumulative error from rounding errors over time 
     1311         ! The following must be bit reproducible over different PE decompositions 
     1312         zantarctica_icesheet_mass_in = glob_sum( antarctica_icesheet_mass_array(:,:) * tmask(:,:,1) ) 
     1313 
     1314         zantarctica_icesheet_mass_in = zantarctica_icesheet_mass_in / zmask_sum 
     1315         antarctica_icesheet_timelapsed = antarctica_icesheet_timelapsed + rdt          
     1316         IF( ABS( zantarctica_icesheet_mass_in - antarctica_icesheet_mass ) > zepsilon ) THEN 
     1317            zantarctica_icesheet_mass_b = antarctica_icesheet_mass 
     1318             
     1319            ! Only update the mass if it has increased 
     1320            IF ( (zantarctica_icesheet_mass_in - antarctica_icesheet_mass) > 0.0 ) THEN 
     1321               antarctica_icesheet_mass = zantarctica_icesheet_mass_in 
     1322            END IF 
     1323             
     1324            IF( zantarctica_icesheet_mass_b /= 0.0 ) & 
     1325          &      antarctica_icesheet_mass_rate_of_change = ( antarctica_icesheet_mass - zantarctica_icesheet_mass_b ) / antarctica_icesheet_timelapsed  
     1326            antarctica_icesheet_timelapsed = 0.0_wp        
     1327         ENDIF 
     1328         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass (kg) read in is ', zantarctica_icesheet_mass_in 
     1329         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass (kg) used is    ', antarctica_icesheet_mass 
     1330         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass rate of change (kg/s) is ', antarctica_icesheet_mass_rate_of_change 
     1331         IF(lwp) WRITE(numout,*) 'Antarctica icesheet seconds lapsed since last change is ', antarctica_icesheet_timelapsed 
     1332      ENDIF 
     1333 
    11121334      ! 
    11131335      CALL wrk_dealloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr ) 
     
    14031625         zsprecip(:,:) = frcv(jpr_snow)%z3(:,:,1)                  ! May need to ensure positive here 
    14041626         ztprecip(:,:) = frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here 
    1405          zemp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:) 
     1627         zemp_tot(:,:) = frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)          
     1628#if defined key_cice 
     1629         IF ( TRIM(sn_rcv_emp%clcat) == 'yes' ) THEN 
     1630            ! zemp_ice is the sum of frcv(jpr_ievp)%z3(:,:,1) over all layers - snow 
     1631            zemp_ice(:,:) = - frcv(jpr_snow)%z3(:,:,1) 
     1632            DO jl=1,jpl 
     1633               zemp_ice(:,:   ) = zemp_ice(:,:) + frcv(jpr_ievp)%z3(:,:,jl) 
     1634            ENDDO 
     1635            ! latent heat coupled for each category in CICE 
     1636            qla_ice(:,:,1:jpl) = - frcv(jpr_ievp)%z3(:,:,1:jpl) * lsub 
     1637         ELSE 
     1638            ! If CICE has multicategories it still expects coupling fields for 
     1639            ! each even if we treat as a single field 
     1640            ! The latent heat flux is split between the ice categories according 
     1641            ! to the fraction of the ice in each category 
     1642            zemp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) 
     1643            WHERE ( zicefr(:,:) /= 0._wp )  
     1644               ztmp(:,:) = 1./zicefr(:,:) 
     1645            ELSEWHERE  
     1646               ztmp(:,:) = 0.e0 
     1647            END WHERE   
     1648            DO jl=1,jpl 
     1649               qla_ice(:,:,jl) = - a_i(:,:,jl) * ztmp(:,:) * frcv(jpr_ievp)%z3(:,:,1) * lsub  
     1650            END DO 
     1651            WHERE ( zicefr(:,:) == 0._wp )  qla_ice(:,:,1) = -frcv(jpr_ievp)%z3(:,:,1) * lsub  
     1652         ENDIF 
     1653#else          
    14061654         zemp_ice(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) 
     1655#endif                   
    14071656            CALL iom_put( 'rain'         , frcv(jpr_rain)%z3(:,:,1)              )   ! liquid precipitation  
    14081657         IF( iom_use('hflx_rain_cea') )   & 
     
    17582007               CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' ) 
    17592008               END SELECT 
     2009            CASE( 'oce and weighted ice' )   ;   ztmp1(:,:) =   tsn(:,:,1,jp_tem) + rt0  
     2010               SELECT CASE( sn_snd_temp%clcat ) 
     2011               CASE( 'yes' )    
     2012           ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl) 
     2013               CASE( 'no' ) 
     2014           ztmp3(:,:,:) = 0.0 
     2015           DO jl=1,jpl 
     2016                     ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl) 
     2017           ENDDO 
     2018               CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' ) 
     2019               END SELECT 
    17602020            CASE( 'mixed oce-ice'        )    
    17612021               ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:)  
     
    17992059         END SELECT 
    18002060         IF( ssnd(jps_fice)%laction )   CALL cpl_snd( jps_fice, isec, ztmp3, info ) 
     2061      ENDIF 
     2062       
     2063      ! Send ice fraction field (first order interpolation), for weighting UM fluxes to be passed to NEMO 
     2064      IF (ssnd(jps_fice1)%laction) THEN 
     2065         SELECT CASE (sn_snd_thick1%clcat) 
     2066         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl) 
     2067         CASE( 'no' )    ;   ztmp3(:,:,1) = fr_i(:,:) 
     2068         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick1%clcat' ) 
     2069    END SELECT 
     2070         CALL cpl_snd (jps_fice1, isec, ztmp3, info) 
    18012071      ENDIF 
    18022072       
     
    18442114         IF( ssnd(jps_hsnw)%laction )   CALL cpl_snd( jps_hsnw, isec, ztmp4, info ) 
    18452115      ENDIF 
     2116      ! 
     2117      ! Send meltpond fields  
     2118      IF( ssnd(jps_a_p)%laction .OR. ssnd(jps_ht_p)%laction ) THEN 
     2119         SELECT CASE( sn_snd_mpnd%cldes)  
     2120         CASE( 'weighted ice' )  
     2121            SELECT CASE( sn_snd_mpnd%clcat )  
     2122            CASE( 'yes' )  
     2123               ztmp3(:,:,1:jpl) =  a_p(:,:,1:jpl) * a_i(:,:,1:jpl)  
     2124               ztmp4(:,:,1:jpl) =  ht_p(:,:,1:jpl) * a_i(:,:,1:jpl)  
     2125            CASE( 'no' )  
     2126               ztmp3(:,:,:) = 0.0  
     2127               ztmp4(:,:,:) = 0.0  
     2128               DO jl=1,jpl  
     2129                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_p(:,:,jpl) * a_i(:,:,jpl)  
     2130                 ztmp4(:,:,1) = ztmp4(:,:,1) + ht_p(:,:,jpl) * a_i(:,:,jpl)  
     2131               ENDDO  
     2132            CASE default    ;   CALL ctl_stop( 'sbc_cpl_mpd: wrong definition of sn_snd_mpnd%clcat' )  
     2133            END SELECT  
     2134         CASE( 'ice only' )     
     2135            ztmp3(:,:,1:jpl) = a_p(:,:,1:jpl)  
     2136            ztmp4(:,:,1:jpl) = ht_p(:,:,1:jpl)  
     2137         END SELECT  
     2138         IF( ssnd(jps_a_p)%laction )   CALL cpl_snd( jps_a_p, isec, ztmp3, info )     
     2139         IF( ssnd(jps_ht_p)%laction )   CALL cpl_snd( jps_ht_p, isec, ztmp4, info )     
     2140         ! 
     2141         ! Send ice effective conductivity 
     2142         SELECT CASE( sn_snd_cond%cldes) 
     2143         CASE( 'weighted ice' )    
     2144            SELECT CASE( sn_snd_cond%clcat ) 
     2145            CASE( 'yes' )    
     2146               ztmp3(:,:,1:jpl) =  kn_ice(:,:,1:jpl) * a_i(:,:,1:jpl) 
     2147            CASE( 'no' ) 
     2148               ztmp3(:,:,:) = 0.0 
     2149               DO jl=1,jpl 
     2150                 ztmp3(:,:,1) = ztmp3(:,:,1) + kn_ice(:,:,jl) * a_i(:,:,jl) 
     2151               ENDDO 
     2152            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_cond%clcat' ) 
     2153            END SELECT 
     2154         CASE( 'ice only' )    
     2155           ztmp3(:,:,1:jpl) = kn_ice(:,:,1:jpl) 
     2156         END SELECT 
     2157         IF( ssnd(jps_kice)%laction )   CALL cpl_snd( jps_kice, isec, ztmp3, info ) 
     2158      ENDIF 
     2159      ! 
    18462160      ! 
    18472161#if defined key_cpl_carbon_cycle 
     
    20232337      IF( ssnd(jps_rnf   )%laction )  CALL cpl_snd( jps_rnf   , isec, RESHAPE ( rnf , (/jpi,jpj,1/) ), info ) 
    20242338      IF( ssnd(jps_taum  )%laction )  CALL cpl_snd( jps_taum  , isec, RESHAPE ( taum, (/jpi,jpj,1/) ), info ) 
    2025  
     2339       
     2340      ztmp1(:,:) = sstfrz(:,:) + rt0 
     2341      IF( ssnd(jps_sstfrz)%laction )  CALL cpl_snd( jps_sstfrz, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info ) 
     2342      ! 
    20262343      CALL wrk_dealloc( jpi,jpj, zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1 ) 
    20272344      CALL wrk_dealloc( jpi,jpj,jpl, ztmp3, ztmp4 ) 
  • branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice_bitcomp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90

    r6362 r6363  
    1515   USE dom_oce         ! ocean space and time domain 
    1616   USE domvvl 
    17    USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic 
     17   USE eosbn2, only : eos_fzp ! Function to calculate freezing point of seawater 
     18   USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic, rt0 
    1819   USE in_out_manager  ! I/O manager 
    1920   USE iom, ONLY : iom_put,iom_use              ! I/O manager library !!Joakim edit 
     
    3738   USE ice_gather_scatter 
    3839   USE ice_calendar, only: dt 
     40# if defined key_cice4 
    3941   USE ice_state, only: aice,aicen,uvel,vvel,vsno,vsnon,vice,vicen 
    40 # if defined key_cice4 
    4142   USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow,  & 
    4243                strocnxT,strocnyT,                               &  
     
    4546                flatn_f,fsurfn_f,fcondtopn_f,                    & 
    4647                uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl,   & 
    47                 swvdr,swvdf,swidr,swidf 
     48                swvdr,swvdf,swidr,swidf,Tf 
    4849   USE ice_therm_vertical, only: calc_Tsfc 
    4950#else 
     51   USE ice_state, only: aice,aicen,uvel,nt_hpnd,trcrn,vvel,vsno,& 
     52                vsnon,vice,vicen,nt_Tsfc 
    5053   USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow,  & 
    5154                strocnxT,strocnyT,                               &  
    52                 sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai,     & 
    53                 fresh_ai,fhocn_ai,fswthru_ai,frzmlt,          & 
     55                sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai,      & 
     56                fresh_ai,fhocn_ai,fswthru_ai,frzmlt,             & 
    5457                flatn_f,fsurfn_f,fcondtopn_f,                    & 
    5558                uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl,   & 
    56                 swvdr,swvdf,swidr,swidf 
    57    USE ice_therm_shared, only: calc_Tsfc 
     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 
    5868#endif 
    5969   USE ice_forcing, only: frcvdr,frcvdf,frcidr,frcidf 
     
    161171      INTEGER, INTENT( in  ) ::   ksbc                ! surface forcing type 
    162172      REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2 
    163       REAL(wp) ::   zcoefu, zcoefv, zcoeff            ! local scalar 
     173      REAL(wp), DIMENSION(:,:,:), POINTER :: ztfrz3d 
    164174      INTEGER  ::   ji, jj, jl, jk                    ! dummy loop indices 
    165175      !!--------------------------------------------------------------------- 
     
    174184      jj_off = INT ( (jpjglo - ny_global) / 2 ) 
    175185 
    176 #if defined key_nemocice_decomp 
    177       ! Pass initial SST from NEMO to CICE so ice is initialised correctly if 
    178       ! there is no restart file. 
    179       ! Values from a CICE restart file would overwrite this 
    180       IF ( .NOT. ln_rstart ) THEN     
    181          CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.)  
    182       ENDIF   
    183 #endif 
    184  
    185 ! Initialize CICE 
     186      ! Initialize CICE 
    186187      CALL CICE_Initialize 
    187188 
    188 ! Do some CICE consistency checks 
     189      ! Do some CICE consistency checks 
    189190      IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN 
    190191         IF ( calc_strair .OR. calc_Tsfc ) THEN 
     
    198199 
    199200 
    200 ! allocate sbc_ice and sbc_cice arrays 
    201       IF( sbc_ice_alloc()      /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate arrays' ) 
     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' ) 
    202203      IF( sbc_ice_cice_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate cice arrays' ) 
    203204 
    204 ! Ensure ocean temperatures are nowhere below freezing if not a NEMO restart 
     205      ! Ensure that no temperature points are below freezing if not a NEMO restart 
    205206      IF( .NOT. ln_rstart ) THEN 
    206          tsn(:,:,:,jp_tem) = MAX (tsn(:,:,:,jp_tem),Tocnfrz) 
     207 
     208         CALL wrk_alloc( jpi,jpj,jpk, ztfrz3d )  
     209         DO jk=1,jpk 
     210            ztfrz3d(:,:,jk) = eos_fzp( tsn(:,:,jk,jp_sal), fsdept_n(:,:,jk) ) 
     211         ENDDO 
     212         tsn(:,:,:,jp_tem) = MAX( tsn(:,:,:,jp_tem), ztfrz3d ) 
    207213         tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) 
    208       ENDIF 
    209  
    210       fr_iu(:,:)=0.0 
    211       fr_iv(:,:)=0.0 
     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.)  
     221#endif 
     222 
     223      ENDIF   
     224 
     225      ! calculate surface freezing temperature and send to CICE 
     226      sstfrz(:,:) = eos_fzp(sss_m(:,:), fsdept_n(:,:,1))  
     227      CALL nemo2cice(sstfrz,Tf, 'T', 1. ) 
    212228 
    213229      CALL cice2nemo(aice,fr_i, 'T', 1. ) 
     
    220236! T point to U point 
    221237! T point to V point 
     238      fr_iu(:,:)=0.0 
     239      fr_iv(:,:)=0.0 
    222240      DO jj=1,jpjm1 
    223241         DO ji=1,jpim1 
     
    343361         CALL nemo2cice(ztmp,stray,'F', -1. ) 
    344362 
     363 
     364! Alex West: From configuration GSI8 onwards, when NEMO is used with CICE in 
     365! HadGEM3 the 'time-travelling ice' coupling approach is used, whereby  
     366! atmosphere-ice fluxes are passed as pseudo-local values, formed by dividing 
     367! gridbox mean fluxes in the UM by future ice concentration obtained through   
     368! OASIS.  This allows for a much more realistic apportionment of energy through 
     369! the ice - and conserves energy. 
     370! Therefore the fluxes are now divided by ice concentration in the coupled 
     371! formulation (jp_purecpl) as well as for jp_flx.  This NEMO branch should only 
     372! be used at UM10.2 onwards (unless an explicit GSI8 UM branch is included), at 
     373! which point the GSI8 UM changes were committed. 
     374 
    345375! Surface downward latent heat flux (CI_5) 
    346376         IF (ksbc == jp_flx) THEN 
     
    348378               ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl) 
    349379            ENDDO 
    350          ELSE 
    351 ! emp_ice is set in sbc_cpl_ice_flx as sublimation-snow 
    352             qla_ice(:,:,1)= - ( emp_ice(:,:)+sprecip(:,:) ) * Lsub 
    353 ! End of temporary code 
    354             DO jj=1,jpj 
    355                DO ji=1,jpi 
    356                   IF (fr_i(ji,jj).eq.0.0) THEN 
    357                      DO jl=1,ncat 
    358                         ztmpn(ji,jj,jl)=0.0 
    359                      ENDDO 
    360                      ! This will then be conserved in CICE 
    361                      ztmpn(ji,jj,1)=qla_ice(ji,jj,1) 
    362                   ELSE 
    363                      DO jl=1,ncat 
    364                         ztmpn(ji,jj,jl)=qla_ice(ji,jj,1)*a_i(ji,jj,jl)/fr_i(ji,jj) 
    365                      ENDDO 
    366                   ENDIF 
    367                ENDDO 
     380         ELSE IF (ksbc == jp_purecpl) THEN 
     381            DO jl=1,ncat 
     382               ztmpn(:,:,jl)=qla_ice(:,:,jl)*a_i(:,:,jl) 
    368383            ENDDO 
     384    ELSE 
     385           !In coupled mode - qla_ice calculated in sbc_cpl for each category 
     386           ztmpn(:,:,1:ncat)=qla_ice(:,:,1:ncat) 
    369387         ENDIF 
     388 
    370389         DO jl=1,ncat 
    371390            CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. ) 
     
    373392! GBM conductive flux through ice (CI_6) 
    374393!  Convert to GBM 
    375             IF (ksbc == jp_flx) THEN 
     394            IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN 
    376395               ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl) 
    377396            ELSE 
     
    382401! GBM surface heat flux (CI_7) 
    383402!  Convert to GBM 
    384             IF (ksbc == jp_flx) THEN 
     403            IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN 
    385404               ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl)  
    386405            ELSE 
     
    442461      CALL nemo2cice(ztmp,frain,'T', 1. )  
    443462 
     463! Recalculate freezing temperature and send to CICE  
     464      sstfrz(:,:)=eos_fzp(sss_m(:,:), fsdept_n(:,:,1))  
     465      CALL nemo2cice(sstfrz,Tf,'T', 1. ) 
     466 
    444467! Freezing/melting potential 
    445468! Calculated over NEMO leapfrog timestep (hence 2*dt) 
    446       nfrzmlt(:,:)=rau0*rcp*fse3t_m(:,:)*(Tocnfrz-sst_m(:,:))/(2.0*dt) 
    447  
    448       ztmp(:,:) = nfrzmlt(:,:) 
    449       CALL nemo2cice(ztmp,frzmlt,'T', 1. ) 
     469      nfrzmlt(:,:)=rau0*rcp*fse3t_m(:,:)*(sstfrz(:,:)-sst_m(:,:))/(2.0*dt)  
     470      CALL nemo2cice(nfrzmlt,frzmlt,'T', 1. ) 
    450471 
    451472! SST  and SSS 
     
    453474      CALL nemo2cice(sst_m,sst,'T', 1. ) 
    454475      CALL nemo2cice(sss_m,sss,'T', 1. ) 
     476 
     477      IF( ksbc == jp_purecpl ) THEN 
     478! Sea ice surface skin temperature 
     479         DO jl=1,ncat 
     480           CALL nemo2cice(tsfc_ice(:,:,jl), trcrn(:,:,nt_tsfc,jl,:),'T',1.) 
     481         ENDDO  
     482      ENDIF 
    455483 
    456484! x comp and y comp of surface ocean current 
     
    730758         CALL cice2nemo(vicen(:,:,jl,:),ht_i(:,:,jl),'T', 1. ) 
    731759      ENDDO 
     760 
     761#if ! defined key_cice4 
     762! Meltpond fraction and depth 
     763      DO jl = 1,ncat 
     764         CALL cice2nemo(apeffn(:,:,jl,:),a_p(:,:,jl),'T', 1. ) 
     765         CALL cice2nemo(trcrn(:,:,nt_hpnd,jl,:),ht_p(:,:,jl),'T', 1. ) 
     766      ENDDO 
     767#endif 
     768 
     769 
     770! If using multilayers thermodynamics in CICE then get top layer temperature 
     771! and effective conductivity        
     772!! When using NEMO with CICE, this change requires use of  
     773!! one of the following two CICE branches: 
     774!! - at CICE5.0,   hadax/r1015_GSI8_with_GSI7 
     775!! - at CICE5.1.2, hadax/vn5.1.2_GSI8 
     776      IF (heat_capacity) THEN 
     777         DO jl = 1,ncat 
     778            CALL cice2nemo(Tn_top(:,:,jl,:),tn_ice(:,:,jl),'T', 1. ) 
     779            CALL cice2nemo(keffn_top(:,:,jl,:),kn_ice(:,:,jl),'T', 1. ) 
     780         ENDDO 
     781! Convert surface temperature to Kelvin 
     782         tn_ice(:,:,:)=tn_ice(:,:,:)+rt0 
     783      ELSE 
     784         tn_ice(:,:,:) = 0.0 
     785         kn_ice(:,:,:) = 0.0 
     786      ENDIF        
     787 
    732788      ! 
    733789      IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_hadgam') 
  • branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice_bitcomp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

    r6362 r6363  
    2626   USE zdfbfr 
    2727   USE fldread         ! read input field at current time step 
    28  
    29  
     28   USE lib_fortran, ONLY: glob_sum 
    3029 
    3130   IMPLICIT NONE 
     
    9089    INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror 
    9190    INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer 
     91    REAL(wp)                     ::   zgreenland_fwfisf_sum, zantarctica_fwfisf_sum 
    9292    REAL(wp)                     ::   rmin 
    9393    REAL(wp)                     ::   zhk 
     
    256256            CALL fld_read ( kt, nn_fsbc, sf_rnfisf   ) 
    257257            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting)  
     258 
     259            IF( lk_oasis) THEN 
     260            ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true 
     261            IF( ln_coupled_iceshelf_fluxes ) THEN 
     262 
     263              ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern 
     264              ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets 
     265              ! to preserve total freshwater conservation in coupled models without an active ice sheet model. 
     266  
     267              ! All related global sums must be done bit reproducibly 
     268               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) ) 
     269 
     270               ! use ABS function because we need to preserve the sign of fwfisf 
     271               WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  & 
     272              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) & 
     273              &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) ) 
     274 
     275               ! check 
     276               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum 
     277 
     278               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) ) 
     279 
     280               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum 
     281 
     282               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) ) 
     283 
     284               ! use ABS function because we need to preserve the sign of fwfisf 
     285               WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) & 
     286              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) & 
     287              &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) ) 
     288       
     289               ! check 
     290               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum 
     291 
     292               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) ) 
     293 
     294               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum 
     295 
     296            ENDIF 
     297            ENDIF 
     298 
    258299            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux 
    259300            stbl(:,:)   = soce 
     
    264305            !CALL fld_read ( kt, nn_fsbc, sf_qisf   ) 
    265306            fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf 
     307 
     308            IF( lk_oasis) THEN 
     309            ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true 
     310            IF( ln_coupled_iceshelf_fluxes ) THEN 
     311 
     312              ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern 
     313              ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets 
     314              ! to preserve total freshwater conservation in coupled models without an active ice sheet model. 
     315 
     316              ! All related global sums must be done bit reproducibly 
     317               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) ) 
     318 
     319               ! use ABS function because we need to preserve the sign of fwfisf 
     320               WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  & 
     321              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) & 
     322              &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) ) 
     323 
     324               ! check 
     325               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum 
     326 
     327               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) ) 
     328 
     329               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum 
     330 
     331               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) ) 
     332 
     333               ! use ABS function because we need to preserve the sign of fwfisf 
     334               WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) & 
     335              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) & 
     336              &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) ) 
     337       
     338               ! check 
     339               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum 
     340 
     341               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) ) 
     342 
     343               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum 
     344 
     345            ENDIF 
     346            ENDIF 
     347 
    266348            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux 
    267349            !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)              ! heat flux 
Note: See TracChangeset for help on using the changeset viewer.