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 6367 – NEMO

Changeset 6367


Ignore:
Timestamp:
2016-03-08T10:08:05+01:00 (8 years ago)
Author:
frrh
Message:

Commit addition of code from
svn merge -r 6340:6343 branches/UKMO/dev_r5518_fix_global_ice

Location:
branches/UKMO/dev_r5518_debug_isf_restart/NEMOGCM
Files:
8 edited

Legend:

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

    r5517 r6367  
    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_debug_isf_restart/NEMOGCM/CONFIG/SHARED/namelist_ref

    r5501 r6367  
    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_debug_isf_restart/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90

    r6366 r6367  
    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_debug_isf_restart/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r6366 r6367  
    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_debug_isf_restart/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r6366 r6367  
    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_debug_isf_restart/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r6366 r6367  
    5252   USE limthd_dh       ! for CALL lim_thd_snwblow 
    5353#endif 
     54   USE lib_fortran, ONLY: glob_sum 
    5455 
    5556   IMPLICIT NONE 
     
    105106   INTEGER, PARAMETER ::   jpr_e3t1st = 41            ! first T level thickness  
    106107   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 
     108   INTEGER, PARAMETER ::   jpr_ts_ice = 43            ! skin temperature of sea-ice (used for melt-ponds) 
     109   INTEGER, PARAMETER ::   jpr_grnm   = 44            ! Greenland ice mass 
     110   INTEGER, PARAMETER ::   jpr_antm   = 45            ! Antarctic ice mass 
     111   INTEGER, PARAMETER ::   jprcv      = 45            ! total number of fields received 
    108112 
    109113   INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere 
     
    149153   ! Received from the atmosphere                     ! 
    150154   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                         
     155   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_ts_ice, sn_rcv_grnm, sn_rcv_antm 
    152156   ! Other namelist parameters                        ! 
    153157   INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data 
     
    216220      REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos 
    217221      !! 
    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 
     222      NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick , sn_snd_crt   , sn_snd_co2,     & 
     223         &                  sn_snd_cond, sn_snd_mpnd  , sn_snd_sstfrz, sn_snd_thick1,                 & 
     224         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,     & 
     225         &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   , sn_rcv_iceflx,  & 
     226         &                  sn_rcv_co2 , sn_rcv_grnm  , sn_rcv_antm  , sn_rcv_ts_ice, nn_cplmodel  ,  & 
     227         &                  ln_usecplmask, ln_coupled_iceshelf_fluxes, rn_greenland_calving_fraction, & 
     228         &                  rn_antarctica_calving_fraction, rn_iceshelf_fluxes_tolerance 
    222229      !!--------------------------------------------------------------------- 
    223230      ! 
     
    258265         WRITE(numout,*)'      runoffs                         = ', TRIM(sn_rcv_rnf%cldes   ), ' (', TRIM(sn_rcv_rnf%clcat   ), ')' 
    259266         WRITE(numout,*)'      calving                         = ', TRIM(sn_rcv_cal%cldes   ), ' (', TRIM(sn_rcv_cal%clcat   ), ')' 
     267         WRITE(numout,*)'      Greenland ice mass              = ', TRIM(sn_rcv_grnm%cldes  ), ' (', TRIM(sn_rcv_grnm%clcat  ), ')' 
     268         WRITE(numout,*)'      Antarctica ice mass             = ', TRIM(sn_rcv_antm%cldes  ), ' (', TRIM(sn_rcv_antm%clcat  ), ')' 
    260269         WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')' 
    261270         WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')' 
     
    271280         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel 
    272281         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask 
     282         WRITE(numout,*)'  ln_coupled_iceshelf_fluxes          = ', ln_coupled_iceshelf_fluxes 
     283         WRITE(numout,*)'  rn_greenland_calving_fraction       = ', rn_greenland_calving_fraction 
     284         WRITE(numout,*)'  rn_antarctica_calving_fraction      = ', rn_antarctica_calving_fraction 
     285         WRITE(numout,*)'  rn_iceshelf_fluxes_tolerance        = ', rn_iceshelf_fluxes_tolerance 
    273286      ENDIF 
    274287 
     
    410423      ! 
    411424      srcv(jpr_cal   )%clname = 'OCalving'   ;   IF( TRIM( sn_rcv_cal%cldes ) == 'coupled' )   srcv(jpr_cal)%laction = .TRUE. 
     425      srcv(jpr_grnm  )%clname = 'OGrnmass'   ;   IF( TRIM( sn_rcv_grnm%cldes ) == 'coupled' )   srcv(jpr_grnm)%laction = .TRUE. 
     426      srcv(jpr_antm  )%clname = 'OAntmass'   ;   IF( TRIM( sn_rcv_antm%cldes ) == 'coupled' )   srcv(jpr_antm)%laction = .TRUE. 
     427 
    412428 
    413429      !                                                      ! ------------------------- ! 
     
    785801      ncpl_qsr_freq = 86400 / ncpl_qsr_freq 
    786802 
     803      IF( ln_coupled_iceshelf_fluxes ) THEN 
     804          ! Crude masks to separate the Antarctic and Greenland icesheets. Obviously something 
     805          ! more complicated could be done if required. 
     806          greenland_icesheet_mask = 0.0 
     807          WHERE( gphit >= 0.0 ) greenland_icesheet_mask = 1.0 
     808          antarctica_icesheet_mask = 0.0 
     809          WHERE( gphit < 0.0 ) antarctica_icesheet_mask = 1.0 
     810 
     811          ! initialise other variables 
     812          greenland_icesheet_mass_array(:,:) = 0.0 
     813          antarctica_icesheet_mass_array(:,:) = 0.0 
     814 
     815          IF( .not. ln_rstart ) THEN 
     816             greenland_icesheet_mass = 0.0  
     817             greenland_icesheet_mass_rate_of_change = 0.0  
     818             greenland_icesheet_timelapsed = 0.0 
     819             antarctica_icesheet_mass = 0.0  
     820             antarctica_icesheet_mass_rate_of_change = 0.0  
     821             antarctica_icesheet_timelapsed = 0.0 
     822          ENDIF 
     823 
     824      ENDIF 
     825 
    787826      CALL wrk_dealloc( jpi,jpj, zacs, zaos ) 
    788827      ! 
     
    846885      INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdttra did not change since nit000) 
    847886      REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars      
     887      REAL(wp) ::   zgreenland_icesheet_mass_in, zantarctica_icesheet_mass_in 
     888      REAL(wp) ::   zgreenland_icesheet_mass_b, zantarctica_icesheet_mass_b 
     889      REAL(wp) ::   zmask_sum, zepsilon       
    848890      REAL(wp) ::   zcoef                  ! temporary scalar 
    849891      REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3 
     
    11101152 
    11111153      ENDIF 
     1154       
     1155      !                                                        ! land ice masses : Greenland 
     1156      zepsilon = rn_iceshelf_fluxes_tolerance 
     1157 
     1158 
     1159      ! See if we need zmask_sum... 
     1160      IF ( srcv(jpr_grnm)%laction .OR. srcv(jpr_antm)%laction ) THEN 
     1161         zmask_sum = glob_sum( tmask(:,:,1) ) 
     1162      ENDIF 
     1163 
     1164      IF( srcv(jpr_grnm)%laction ) THEN 
     1165         greenland_icesheet_mass_array(:,:) = frcv(jpr_grnm)%z3(:,:,1) 
     1166         ! take average over ocean points of input array to avoid cumulative error over time 
     1167 
     1168         ! The following must be bit reproducible over different PE decompositions 
     1169         zgreenland_icesheet_mass_in = glob_sum( greenland_icesheet_mass_array(:,:) * tmask(:,:,1) ) 
     1170 
     1171         zgreenland_icesheet_mass_in = zgreenland_icesheet_mass_in / zmask_sum 
     1172         greenland_icesheet_timelapsed = greenland_icesheet_timelapsed + rdt          
     1173         IF( ABS( zgreenland_icesheet_mass_in - greenland_icesheet_mass ) > zepsilon ) THEN 
     1174            zgreenland_icesheet_mass_b = greenland_icesheet_mass 
     1175             
     1176            ! Only update the mass if it has increased 
     1177            IF ( (zgreenland_icesheet_mass_in - greenland_icesheet_mass) > 0.0 ) THEN 
     1178               greenland_icesheet_mass = zgreenland_icesheet_mass_in 
     1179            ENDIF 
     1180             
     1181            IF( zgreenland_icesheet_mass_b /= 0.0 ) & 
     1182           &     greenland_icesheet_mass_rate_of_change = ( greenland_icesheet_mass - zgreenland_icesheet_mass_b ) / greenland_icesheet_timelapsed  
     1183            greenland_icesheet_timelapsed = 0.0_wp        
     1184         ENDIF 
     1185         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass (kg) read in is ', zgreenland_icesheet_mass_in 
     1186         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass (kg) used is    ', greenland_icesheet_mass 
     1187         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass rate of change (kg/s) is ', greenland_icesheet_mass_rate_of_change 
     1188         IF(lwp) WRITE(numout,*) 'Greenland icesheet seconds lapsed since last change is ', greenland_icesheet_timelapsed 
     1189      ENDIF 
     1190 
     1191      !                                                        ! land ice masses : Antarctica 
     1192      IF( srcv(jpr_antm)%laction ) THEN 
     1193         antarctica_icesheet_mass_array(:,:) = frcv(jpr_antm)%z3(:,:,1) 
     1194         ! take average over ocean points of input array to avoid cumulative error from rounding errors over time 
     1195         ! The following must be bit reproducible over different PE decompositions 
     1196         zantarctica_icesheet_mass_in = glob_sum( antarctica_icesheet_mass_array(:,:) * tmask(:,:,1) ) 
     1197 
     1198         zantarctica_icesheet_mass_in = zantarctica_icesheet_mass_in / zmask_sum 
     1199         antarctica_icesheet_timelapsed = antarctica_icesheet_timelapsed + rdt          
     1200         IF( ABS( zantarctica_icesheet_mass_in - antarctica_icesheet_mass ) > zepsilon ) THEN 
     1201            zantarctica_icesheet_mass_b = antarctica_icesheet_mass 
     1202             
     1203            ! Only update the mass if it has increased 
     1204            IF ( (zantarctica_icesheet_mass_in - antarctica_icesheet_mass) > 0.0 ) THEN 
     1205               antarctica_icesheet_mass = zantarctica_icesheet_mass_in 
     1206            END IF 
     1207             
     1208            IF( zantarctica_icesheet_mass_b /= 0.0 ) & 
     1209          &      antarctica_icesheet_mass_rate_of_change = ( antarctica_icesheet_mass - zantarctica_icesheet_mass_b ) / antarctica_icesheet_timelapsed  
     1210            antarctica_icesheet_timelapsed = 0.0_wp        
     1211         ENDIF 
     1212         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass (kg) read in is ', zantarctica_icesheet_mass_in 
     1213         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass (kg) used is    ', antarctica_icesheet_mass 
     1214         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass rate of change (kg/s) is ', antarctica_icesheet_mass_rate_of_change 
     1215         IF(lwp) WRITE(numout,*) 'Antarctica icesheet seconds lapsed since last change is ', antarctica_icesheet_timelapsed 
     1216      ENDIF 
     1217 
    11121218      ! 
    11131219      CALL wrk_dealloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr ) 
  • branches/UKMO/dev_r5518_debug_isf_restart/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90

    r6366 r6367  
    161161      INTEGER, INTENT( in  ) ::   ksbc                ! surface forcing type 
    162162      REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2 
    163       REAL(wp) ::   zcoefu, zcoefv, zcoeff            ! local scalar 
     163      REAL(wp), DIMENSION(:,:,:), POINTER :: ztfrz3d 
    164164      INTEGER  ::   ji, jj, jl, jk                    ! dummy loop indices 
    165165      !!--------------------------------------------------------------------- 
     
    174174      jj_off = INT ( (jpjglo - ny_global) / 2 ) 
    175175 
    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 
     176      ! Initialize CICE 
    186177      CALL CICE_Initialize 
    187178 
    188 ! Do some CICE consistency checks 
     179      ! Do some CICE consistency checks 
    189180      IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN 
    190181         IF ( calc_strair .OR. calc_Tsfc ) THEN 
     
    198189 
    199190 
    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' ) 
     191      ! allocate sbc_ice and sbc_cice arrays 
     192      IF( sbc_ice_alloc()      /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_alloc : unable to allocate arrays' ) 
    202193      IF( sbc_ice_cice_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate cice arrays' ) 
    203194 
    204 ! Ensure ocean temperatures are nowhere below freezing if not a NEMO restart 
     195      ! Ensure that no temperature points are below freezing if not a NEMO restart 
    205196      IF( .NOT. ln_rstart ) THEN 
    206          tsn(:,:,:,jp_tem) = MAX (tsn(:,:,:,jp_tem),Tocnfrz) 
     197 
     198         CALL wrk_alloc( jpi,jpj,jpk, ztfrz3d )  
     199         DO jk=1,jpk 
     200            ztfrz3d(:,:,jk) = eos_fzp( tsn(:,:,jk,jp_sal), fsdept_n(:,:,jk) ) 
     201         ENDDO 
     202         tsn(:,:,:,jp_tem) = MAX( tsn(:,:,:,jp_tem), ztfrz3d ) 
    207203         tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) 
    208       ENDIF 
    209  
    210       fr_iu(:,:)=0.0 
    211       fr_iv(:,:)=0.0 
     204         CALL wrk_dealloc( jpi,jpj,jpk, ztfrz3d )  
     205 
     206#if defined key_nemocice_decomp 
     207         ! Pass initial SST from NEMO to CICE so ice is initialised correctly if 
     208         ! there is no restart file. 
     209         ! Values from a CICE restart file would overwrite this 
     210         CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.)  
     211#endif 
     212 
     213      ENDIF   
     214 
     215      ! calculate surface freezing temperature and send to CICE 
     216      sstfrz(:,:) = eos_fzp(sss_m(:,:), fsdept_n(:,:,1))  
     217      CALL nemo2cice(sstfrz,Tf, 'T', 1. ) 
    212218 
    213219      CALL cice2nemo(aice,fr_i, 'T', 1. ) 
     
    220226! T point to U point 
    221227! T point to V point 
     228      fr_iu(:,:)=0.0 
     229      fr_iv(:,:)=0.0 
    222230      DO jj=1,jpjm1 
    223231         DO ji=1,jpim1 
     
    348356               ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl) 
    349357            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 
     358         ELSE IF (ksbc == jp_purecpl) THEN 
     359            DO jl=1,ncat 
     360               ztmpn(:,:,jl)=qla_ice(:,:,jl)*a_i(:,:,jl) 
    368361            ENDDO 
     362    ELSE 
     363           !In coupled mode - qla_ice calculated in sbc_cpl for each category 
     364           ztmpn(:,:,1:ncat)=qla_ice(:,:,1:ncat) 
    369365         ENDIF 
    370366         DO jl=1,ncat 
     
    454450      CALL nemo2cice(sss_m,sss,'T', 1. ) 
    455451 
     452      IF( ksbc == jp_purecpl ) THEN 
     453! Sea ice surface skin temperature 
     454         DO jl=1,ncat 
     455           CALL nemo2cice(tsfc_ice(:,:,jl), trcrn(:,:,nt_tsfc,jl,:),'T',1.) 
     456         ENDDO  
     457      ENDIF 
     458 
    456459! x comp and y comp of surface ocean current 
    457460! U point to F point 
  • branches/UKMO/dev_r5518_debug_isf_restart/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

    r6366 r6367  
    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.