Changeset 6200

Timestamp:

2016-01-04T11:49:32+01:00 (8 years ago)

Author:

frrh

Message:

Merge in branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice@5797
and MY medusa interface, resolve conflicts in sbccpl and, mysteriously
in sbcice_cice.F90 which frankly should not occur since I am doing nothing in
here!

Location:

branches/UKMO/dev_r5518_med_test/NEMOGCM

Files:

• CONFIG/SHARED/field_def.xml (modified) (1 diff)
• CONFIG/SHARED/namelist_ref (modified) (1 diff)
• NEMO/OPA_SRC/ICB/icbclv.F90 (modified) (3 diffs)
• NEMO/OPA_SRC/IOM/restart.F90 (modified) (3 diffs)
• NEMO/OPA_SRC/SBC/sbc_oce.F90 (modified) (3 diffs)
• NEMO/OPA_SRC/SBC/sbccpl.F90 (modified) (18 diffs)
• NEMO/OPA_SRC/SBC/sbcice_cice.F90 (modified) (5 diffs)
• NEMO/OPA_SRC/SBC/sbcisf.F90 (modified) (4 diffs)

branches/UKMO/dev_r5518_med_test/NEMOGCM/CONFIG/SHARED/field_def.xml

@@ -193 +193 @@
 
          <!-- * variable related to ice shelf forcing * -->
+         <field id="berg_calve"   long_name="Iceberg calving"                               unit="kg/m2/s"  />
          <field id="fwfisf"       long_name="Ice shelf melting"                             unit="kg/m2/s"  />
          <field id="qisf"         long_name="Ice Shelf Heat Flux"                           unit="W/m2"     />

branches/UKMO/dev_r5518_med_test/NEMOGCM/CONFIG/SHARED/namelist_ref

@@ -378 +378 @@
    ln_usecplmask = .false.   !  use a coupling mask file to merge data received from several models
                              !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel)
+   ln_coupled_iceshelf_fluxes = .false. ! If true use rate of change of mass of Greenland and Antarctic icesheets to set the 
+                                        ! combined magnitude of the iceberg calving and iceshelf melting freshwater fluxes.
+   rn_greenland_calving_fraction = 0.5  ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
+   rn_antarctica_calving_fraction = 0.5 ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
+   rn_iceshelf_fluxes_tolerance = 1e-6  ! Fractional threshold for detecting differences in icesheet masses (must be positive definite).
 /
 !-----------------------------------------------------------------------

branches/UKMO/dev_r5518_med_test/NEMOGCM/NEMO/OPA_SRC/ICB/icbclv.F90

@@ -25 +25 @@
    USE icbutl         ! iceberg utility routines
 
+   USE sbc_oce        ! for icesheet freshwater input variables
+   USE in_out_manager
+   USE iom
+
    IMPLICIT NONE
    PRIVATE
@@ -48 +52 @@
       !
       REAL(wp)                        :: zcalving_used, zdist, zfact
+      REAL(wp)                        :: zgreenland_calving_sum, zantarctica_calving_sum
       INTEGER                         :: jn, ji, jj                    ! loop counters
       INTEGER                         :: imx                           ! temporary integer for max berg class
@@ -59 +64 @@
       zfact = ( (1000._wp)**3 / ( NINT(rday) * nyear_len(1) ) ) * 850._wp
       berg_grid%calving(:,:) = src_calving(:,:) * tmask_i(:,:) * zfact
+
+      IF( lk_oasis) THEN
+      ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+      IF( ln_coupled_iceshelf_fluxes ) THEN
+
+        ! Adjust total calving rates so that sum of iceberg calving and iceshelf melting in the northern
+        ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets
+        ! to preserve total freshwater conservation in coupled models without an active ice sheet model.
+
+         zgreenland_calving_sum = SUM( berg_grid%calving(:,:) * greenland_icesheet_mask(:,:) )
+         IF( lk_mpp ) CALL mpp_sum( zgreenland_calving_sum )
+         WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                                 &
+        &    berg_grid%calving(:,:) = berg_grid%calving(:,:) * greenland_icesheet_mass_rate_of_change * rn_greenland_calving_fraction &
+        &                                     / ( zgreenland_calving_sum + 1.0e-10_wp )
+
+         ! check
+         IF(lwp) WRITE(numout, *) 'Greenland iceberg calving climatology (kg/s) : ',zgreenland_calving_sum
+         zgreenland_calving_sum = SUM( berg_grid%calving(:,:) * greenland_icesheet_mask(:,:) )
+         IF( lk_mpp ) CALL mpp_sum( zgreenland_calving_sum )
+         IF(lwp) WRITE(numout, *) 'Greenland iceberg calving adjusted value (kg/s) : ',zgreenland_calving_sum
+
+         zantarctica_calving_sum = SUM( berg_grid%calving(:,:) * antarctica_icesheet_mask(:,:) )
+         IF( lk_mpp ) CALL mpp_sum( zantarctica_calving_sum )
+         WHERE( antarctica_icesheet_mask(:,:) == 1.0 )                                                                              &
+         berg_grid%calving(:,:) = berg_grid%calving(:,:) * antarctica_icesheet_mass_rate_of_change * rn_antarctica_calving_fraction &
+        &                           / ( zantarctica_calving_sum + 1.0e-10_wp )
+ 
+         ! check
+         IF(lwp) WRITE(numout, *) 'Antarctica iceberg calving climatology (kg/s) : ',zantarctica_calving_sum
+         zantarctica_calving_sum = SUM( berg_grid%calving(:,:) * antarctica_icesheet_mask(:,:) )
+         IF( lk_mpp ) CALL mpp_sum( zantarctica_calving_sum )
+         IF(lwp) WRITE(numout, *) 'Antarctica iceberg calving adjusted value (kg/s) : ',zantarctica_calving_sum
+
+      ENDIF
+      ENDIF
+   
+      CALL iom_put( 'berg_calve', berg_grid%calving(:,:) )
 
       ! Heat in units of W/m2, and mask (just in case)

branches/UKMO/dev_r5518_med_test/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

@@ -24 +24 @@
    USE trdmxl_oce      ! ocean active mixed layer tracers trends variables
    USE divcur          ! hor. divergence and curl      (div & cur routines)
+   USE sbc_oce         ! for icesheet freshwater input variables
 
    IMPLICIT NONE
@@ -145 +146 @@
                      CALL iom_rstput( kt, nitrst, numrow, 'rhd'    , rhd       )
 #endif
+                     IF( lk_oasis) THEN
+                     ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+                     IF( ln_coupled_iceshelf_fluxes ) THEN
+                        CALL iom_rstput( kt, nitrst, numrow, 'greenland_icesheet_mass', greenland_icesheet_mass )
+                        CALL iom_rstput( kt, nitrst, numrow, 'greenland_icesheet_timelapsed', greenland_icesheet_timelapsed )
+                        CALL iom_rstput( kt, nitrst, numrow, 'greenland_icesheet_mass_roc', greenland_icesheet_mass_rate_of_change )
+                        CALL iom_rstput( kt, nitrst, numrow, 'antarctica_icesheet_mass', antarctica_icesheet_mass )
+                        CALL iom_rstput( kt, nitrst, numrow, 'antarctica_icesheet_timelapsed', antarctica_icesheet_timelapsed )
+                        CALL iom_rstput( kt, nitrst, numrow, 'antarctica_icesheet_mass_roc', antarctica_icesheet_mass_rate_of_change )
+                     ENDIF
+                     ENDIF
+
       IF( kt == nitrst ) THEN
          CALL iom_close( numrow )     ! close the restart file (only at last time step)
@@ -258 +271 @@
 #endif
       !
+      IF( iom_varid( numror, 'greenland_icesheet_mass', ldstop = .FALSE. ) > 0 )   THEN
+         CALL iom_get( numror, 'greenland_icesheet_mass', greenland_icesheet_mass )
+         CALL iom_get( numror, 'greenland_icesheet_timelapsed', greenland_icesheet_timelapsed )
+         CALL iom_get( numror, 'greenland_icesheet_mass_roc', greenland_icesheet_mass_rate_of_change )
+      ELSE
+         greenland_icesheet_mass = 0.0 
+         greenland_icesheet_mass_rate_of_change = 0.0 
+         greenland_icesheet_timelapsed = 0.0
+      ENDIF
+      IF( iom_varid( numror, 'antarctica_icesheet_mass', ldstop = .FALSE. ) > 0 )   THEN
+         CALL iom_get( numror, 'antarctica_icesheet_mass', antarctica_icesheet_mass )
+         CALL iom_get( numror, 'antarctica_icesheet_timelapsed', antarctica_icesheet_timelapsed )
+         CALL iom_get( numror, 'antarctica_icesheet_mass_roc', antarctica_icesheet_mass_rate_of_change )
+      ELSE
+         antarctica_icesheet_mass = 0.0 
+         antarctica_icesheet_mass_rate_of_change = 0.0 
+         antarctica_icesheet_timelapsed = 0.0
+      ENDIF
+
       IF( neuler == 0 ) THEN                                  ! Euler restart (neuler=0)
          tsb  (:,:,:,:) = tsn  (:,:,:,:)                             ! all before fields set to now values

branches/UKMO/dev_r5518_med_test/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

@@ -125 +125 @@
 #endif
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   greenland_icesheet_mass_array, greenland_icesheet_mask
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   antarctica_icesheet_mass_array, antarctica_icesheet_mask
 
    !!----------------------------------------------------------------------
@@ -137 +139 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3t_m     !: mean (nn_fsbc time-step) sea surface layer thickness       [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frq_m     !: mean (nn_fsbc time-step) fraction of solar net radiation absorbed in the 1st T level [-]
+   
+   !!----------------------------------------------------------------------
+   !!  Surface scalars of total ice sheet mass for Greenland and Antarctica, 
+   !! passed from atmosphere to be converted to dvol and hence a freshwater 
+   !! flux  by using old values. New values are saved in the dump, to become
+   !! old values next coupling timestep. Freshwater fluxes split between 
+   !! sub iceshelf melting and iceberg calving, scalled to flux per second
+   !!----------------------------------------------------------------------
+   
+   REAL(wp), PUBLIC  :: greenland_icesheet_mass, greenland_icesheet_mass_rate_of_change, greenland_icesheet_timelapsed 
+   REAL(wp), PUBLIC  :: antarctica_icesheet_mass, antarctica_icesheet_mass_rate_of_change, antarctica_icesheet_timelapsed
+
+   ! sbccpl namelist parameters associated with icesheet freshwater input code. Included here rather than in sbccpl.F90 to 
+   ! avoid circular dependencies.
+   LOGICAL, PUBLIC     ::   ln_coupled_iceshelf_fluxes     ! If true use rate of change of mass of Greenland and Antarctic icesheets to set the 
+                                                           ! combined magnitude of the iceberg calving and iceshelf melting freshwater fluxes.
+   REAL(wp), PUBLIC    ::   rn_greenland_calving_fraction  ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
+   REAL(wp), PUBLIC    ::   rn_antarctica_calving_fraction ! Set fraction of total freshwater flux for iceberg calving - remainder goes to iceshelf melting.
+   REAL(wp), PUBLIC    ::   rn_iceshelf_fluxes_tolerance   ! Absolute tolerance for detecting differences in icesheet masses. 
 
    !! * Substitutions
@@ -172 +193 @@
          &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) ,      &
          &      ssv_m  (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) )
+      ALLOCATE( greenland_icesheet_mass_array(jpi,jpj) , antarctica_icesheet_mass_array(jpi,jpj) )
+      ALLOCATE( greenland_icesheet_mask(jpi,jpj) , antarctica_icesheet_mask(jpi,jpj) )
          !
 #if defined key_vvl

branches/UKMO/dev_r5518_med_test/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

@@ -52 +52 @@
    USE limthd_dh       ! for CALL lim_thd_snwblow
 #endif
+#if defined key_medusa
+! RSRH Temporarily commented out until we get MEDUSA code   USE trc
+#endif
+
+
+
 
    IMPLICIT NONE
@@ -105 +111 @@
    INTEGER, PARAMETER ::   jpr_e3t1st = 41            ! first T level thickness 
    INTEGER, PARAMETER ::   jpr_fraqsr = 42            ! fraction of solar net radiation absorbed in the first ocean level
-   INTEGER, PARAMETER ::   jprcv      = 42            ! total number of fields received
+   INTEGER, PARAMETER ::   jpr_ts_ice = 43            ! skin temperature of sea-ice (used for melt-ponds)
+   INTEGER, PARAMETER ::   jpr_grnm   = 44            ! Greenland ice mass
+   INTEGER, PARAMETER ::   jpr_antm   = 45            ! Antarctic ice mass
+   INTEGER, PARAMETER ::   jpr_atm_pco2 = 46          ! Incoming atm CO2 flux
+   INTEGER, PARAMETER ::   jpr_atm_dust = 47          ! Incoming atm aggregate dust 
+   INTEGER, PARAMETER ::   jprcv      = 47            ! total number of fields received
 
    INTEGER, PARAMETER ::   jps_fice   =  1            ! ice fraction sent to the atmosphere
@@ -135 +146 @@
    INTEGER, PARAMETER ::   jps_e3t1st = 27            ! first level depth (vvl)
    INTEGER, PARAMETER ::   jps_fraqsr = 28            ! fraction of solar net radiation absorbed in the first ocean level
-   INTEGER, PARAMETER ::   jpsnd      = 28            ! total number of fields sended
+   INTEGER, PARAMETER ::   jps_bio_co2 = 29           ! MEDUSA air-sea CO2 flux in
+   INTEGER, PARAMETER ::   jps_bio_dms = 30           ! MEDUSA DMS surface concentration in
+   INTEGER, PARAMETER ::   jpsnd      = 30            ! total number of fields sent
+
+   REAL(wp), PARAMETER :: dms_unit_conv = 1.0e+6      ! Coversion factor to get outgong DMS in standard units for coupling
+                                                 ! i.e. specifically nmol/L (= umol/m3)
 
    !                                                         !!** namelist namsbc_cpl **
@@ -146 +162 @@
    END TYPE FLD_C
    ! Send to the atmosphere                           !
-   TYPE(FLD_C) ::   sn_snd_temp, sn_snd_alb, sn_snd_thick, sn_snd_crt, sn_snd_co2                        
+   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
+   TYPE(FLD_C) ::   sn_snd_bio_co2, sn_snd_bio_dms                       
    ! Received from the atmosphere                     !
    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
-   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2                        
+   TYPE(FLD_C) ::   sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_ts_ice, sn_rcv_grnm, sn_rcv_antm
+   TYPE(FLD_C) ::   sn_rcv_atm_pco2, sn_rcv_atm_dust                         
    ! Other namelist parameters                        !
    INTEGER     ::   nn_cplmodel            ! Maximum number of models to/from which NEMO is potentialy sending/receiving data
@@ -163 +181 @@
 
    INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument
+
 
    !! Substitution
@@ -216 +235 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos
       !!
-      NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2,      &
-         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr,      &
-         &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf  , sn_rcv_cal   , sn_rcv_iceflx,   &
-         &                  sn_rcv_co2 , nn_cplmodel  , ln_usecplmask
+      NAMELIST/namsbc_cpl/  sn_snd_temp, sn_snd_alb   , sn_snd_thick , sn_snd_crt   , sn_snd_co2,     &
+         &                  sn_snd_bio_co2, sn_snd_bio_dms,                                           &
+         &                  sn_snd_cond, sn_snd_mpnd  , sn_snd_sstfrz, sn_snd_thick1,                 &
+         &                  sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau   , sn_rcv_dqnsdt, sn_rcv_qsr,     &
+         &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   , sn_rcv_iceflx,  &
+         &                  sn_rcv_atm_pco2, sn_rcv_atm_dust,                                         &
+         &                  sn_rcv_co2 , sn_rcv_grnm  , sn_rcv_antm  , sn_rcv_ts_ice, nn_cplmodel  ,  &
+         &                  ln_usecplmask, ln_coupled_iceshelf_fluxes, rn_greenland_calving_fraction, &
+         &                  rn_antarctica_calving_fraction, rn_iceshelf_fluxes_tolerance
+
       !!---------------------------------------------------------------------
       !
@@ -245 +270 @@
       ENDIF
       IF( lwp .AND. ln_cpl ) THEN                        ! control print
-         WRITE(numout,*)'  received fields (mutiple ice categogies)'
+         WRITE(numout,*)'  received fields (mutiple ice categories)'
          WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')'
          WRITE(numout,*)'      stress module                   = ', TRIM(sn_rcv_taumod%cldes), ' (', TRIM(sn_rcv_taumod%clcat), ')'
@@ -260 +285 @@
          WRITE(numout,*)'      sea ice heat fluxes             = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')'
          WRITE(numout,*)'      atm co2                         = ', TRIM(sn_rcv_co2%cldes   ), ' (', TRIM(sn_rcv_co2%clcat   ), ')'
+         WRITE(numout,*)'      atm pco2                        = ', TRIM(sn_rcv_atm_pco2%cldes), ' (', TRIM(sn_rcv_atm_pco2%clcat), ')'
+         WRITE(numout,*)'      atm dust                        = ', TRIM(sn_rcv_atm_dust%cldes), ' (', TRIM(sn_rcv_atm_dust%clcat), ')'
          WRITE(numout,*)'  sent fields (multiple ice categories)'
          WRITE(numout,*)'      surface temperature             = ', TRIM(sn_snd_temp%cldes  ), ' (', TRIM(sn_snd_temp%clcat  ), ')'
@@ -269 +296 @@
          WRITE(numout,*)'                      - mesh          = ', sn_snd_crt%clvgrd
          WRITE(numout,*)'      oce co2 flux                    = ', TRIM(sn_snd_co2%cldes   ), ' (', TRIM(sn_snd_co2%clcat   ), ')'
+         WRITE(numout,*)'      bio co2 flux                    = ', TRIM(sn_snd_bio_co2%cldes), ' (', TRIM(sn_snd_bio_co2%clcat), ')'
+         WRITE(numout,*)'      bio dms flux                    = ', TRIM(sn_snd_bio_dms%cldes), ' (', TRIM(sn_snd_bio_dms%clcat), ')'
          WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel
          WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask
+         WRITE(numout,*)'  ln_coupled_iceshelf_fluxes          = ', ln_coupled_iceshelf_fluxes
+         WRITE(numout,*)'  rn_greenland_calving_fraction       = ', rn_greenland_calving_fraction
+         WRITE(numout,*)'  rn_antarctica_calving_fraction      = ', rn_antarctica_calving_fraction
+         WRITE(numout,*)'  rn_iceshelf_fluxes_tolerance        = ', rn_iceshelf_fluxes_tolerance
       ENDIF
 
@@ -474 +507 @@
       !                                                      ! ------------------------- !
       srcv(jpr_co2 )%clname = 'O_AtmCO2'   ;   IF( TRIM(sn_rcv_co2%cldes   ) == 'coupled' )    srcv(jpr_co2 )%laction = .TRUE.
+
+
+      !                                                      ! --------------------------------------- !    
+      !                                                      ! Incoming CO2 and DUST fluxes for MEDUSA !
+      !                                                      ! --------------------------------------- !  
+      srcv(jpr_atm_pco2)%clname = 'OATMPCO2'
+
+      IF (TRIM(sn_rcv_atm_pco2%cldes) == 'coupled') THEN
+        srcv(jpr_atm_pco2)%laction = .TRUE.
+      END IF
+               
+      srcv(jpr_atm_dust)%clname = 'OATMDUST'   
+      IF (TRIM(sn_rcv_atm_dust%cldes) == 'coupled')  THEN
+        srcv(jpr_atm_dust)%laction = .TRUE.
+      END IF
+    
       !                                                      ! ------------------------- !
       !                                                      !   topmelt and botmelt     !   
@@ -789 +838 @@
       ncpl_qsr_freq = 86400 / ncpl_qsr_freq
 
+      IF( ln_coupled_iceshelf_fluxes ) THEN
+          ! Crude masks to separate the Antarctic and Greenland icesheets. Obviously something
+          ! more complicated could be done if required.
+          greenland_icesheet_mask = 0.0
+          WHERE( gphit >= 0.0 ) greenland_icesheet_mask = 1.0
+          antarctica_icesheet_mask = 0.0
+          WHERE( gphit < 0.0 ) antarctica_icesheet_mask = 1.0
+
+          ! initialise other variables
+          greenland_icesheet_mass_array(:,:) = 0.0
+          antarctica_icesheet_mass_array(:,:) = 0.0
+
+          IF( .not. ln_rstart ) THEN
+             greenland_icesheet_mass = 0.0 
+             greenland_icesheet_mass_rate_of_change = 0.0 
+             greenland_icesheet_timelapsed = 0.0
+             antarctica_icesheet_mass = 0.0 
+             antarctica_icesheet_mass_rate_of_change = 0.0 
+             antarctica_icesheet_timelapsed = 0.0
+          ENDIF
+
+      ENDIF
+
       CALL wrk_dealloc( jpi,jpj, zacs, zaos )
       !
@@ -851 +923 @@
       INTEGER  ::   ikchoix
       REAL(wp) ::   zcumulneg, zcumulpos   ! temporary scalars     
+      REAL(wp) ::   zgreenland_icesheet_mass_in, zantarctica_icesheet_mass_in
+      REAL(wp) ::   zgreenland_icesheet_mass_b, zantarctica_icesheet_mass_b
+      REAL(wp) ::   zmask_sum, zepsilon      
       REAL(wp) ::   zcoef                  ! temporary scalar
       REAL(wp) ::   zrhoa  = 1.22          ! Air density kg/m3
@@ -857 +932 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   ztx, zty, zmsk, zemp, zqns, zqsr, ztx2, zty2
       !!----------------------------------------------------------------------
+
+      ! RSRH temporary arrays for testing, just to recieve incoming MEDUSA related fields
+      ! until we know where they need to go.
+      REAL(wp), ALLOCATABLE :: atm_pco2(:,:)
+      REAL(wp), ALLOCATABLE :: atm_dust(:,:)
+
       !
       IF( nn_timing == 1 )  CALL timing_start('sbc_cpl_rcv')
@@ -1024 +1105 @@
 #endif
 
+#if defined key_medusa
+      ! RSRH Allocate temporary arrays to receive incoming fields during testing
+      ALLOCATE(atm_pco2(jpi,jpj))
+      ALLOCATE(atm_dust(jpi,jpj))
+
+      IF( srcv(jpr_atm_pco2)%laction) atm_pco2(:,:) = frcv(jpr_atm_pco2)%z3(:,:,1)
+      IF( srcv(jpr_atm_dust)%laction) atm_dust(:,:) = frcv(jpr_atm_dust)%z3(:,:,1)
+     
+      ! RSRH Deallocate temporary arrays.
+      DEALLOCATE(atm_pco2)
+      DEALLOCATE(atm_dust)
+#endif
+
+
+
+
       !  Fields received by SAS when OASIS coupling
       !  (arrays no more filled at sbcssm stage)
@@ -1137 +1234 @@
 
       ENDIF
+      
+      !                                                        ! land ice masses : Greenland
+      zepsilon = rn_iceshelf_fluxes_tolerance
+
+      IF( srcv(jpr_grnm)%laction ) THEN
+         greenland_icesheet_mass_array(:,:) = frcv(jpr_grnm)%z3(:,:,1)
+         ! take average over ocean points of input array to avoid cumulative error over time
+         zgreenland_icesheet_mass_in = SUM( greenland_icesheet_mass_array(:,:) * tmask(:,:,1) )
+         IF(lk_mpp) CALL mpp_sum( zgreenland_icesheet_mass_in )
+         zmask_sum = SUM( tmask(:,:,1) )
+         IF(lk_mpp) CALL mpp_sum( zmask_sum ) 
+         zgreenland_icesheet_mass_in = zgreenland_icesheet_mass_in / zmask_sum
+         greenland_icesheet_timelapsed = greenland_icesheet_timelapsed + rdt         
+         IF( ABS( zgreenland_icesheet_mass_in - greenland_icesheet_mass ) > zepsilon ) THEN
+            zgreenland_icesheet_mass_b = greenland_icesheet_mass
+            
+            ! Only update the mass if it has increased
+            IF ( (zgreenland_icesheet_mass_in - greenland_icesheet_mass) > 0.0 ) THEN
+               greenland_icesheet_mass = zgreenland_icesheet_mass_in
+            ENDIF
+            
+            IF( zgreenland_icesheet_mass_b /= 0.0 ) &
+           &     greenland_icesheet_mass_rate_of_change = ( greenland_icesheet_mass - zgreenland_icesheet_mass_b ) / greenland_icesheet_timelapsed 
+            greenland_icesheet_timelapsed = 0.0_wp       
+         ENDIF
+         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass (kg) read in is ', zgreenland_icesheet_mass_in
+         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass (kg) used is    ', greenland_icesheet_mass
+         IF(lwp) WRITE(numout,*) 'Greenland icesheet mass rate of change (kg/s) is ', greenland_icesheet_mass_rate_of_change
+         IF(lwp) WRITE(numout,*) 'Greenland icesheet seconds lapsed since last change is ', greenland_icesheet_timelapsed
+      ENDIF
+
+      !                                                        ! land ice masses : Antarctica
+      IF( srcv(jpr_antm)%laction ) THEN
+         antarctica_icesheet_mass_array(:,:) = frcv(jpr_antm)%z3(:,:,1)
+         ! take average over ocean points of input array to avoid cumulative error from rounding errors over time
+         zantarctica_icesheet_mass_in = SUM( antarctica_icesheet_mass_array(:,:) * tmask(:,:,1) )
+         IF(lk_mpp) CALL mpp_sum( zantarctica_icesheet_mass_in )
+         zmask_sum = SUM( tmask(:,:,1) )
+         IF(lk_mpp) CALL mpp_sum( zmask_sum ) 
+         zantarctica_icesheet_mass_in = zantarctica_icesheet_mass_in / zmask_sum
+         antarctica_icesheet_timelapsed = antarctica_icesheet_timelapsed + rdt         
+         IF( ABS( zantarctica_icesheet_mass_in - antarctica_icesheet_mass ) > zepsilon ) THEN
+            zantarctica_icesheet_mass_b = antarctica_icesheet_mass
+            
+            ! Only update the mass if it has increased
+            IF ( (zantarctica_icesheet_mass_in - antarctica_icesheet_mass) > 0.0 ) THEN
+               antarctica_icesheet_mass = zantarctica_icesheet_mass_in
+            END IF
+            
+            IF( zantarctica_icesheet_mass_b /= 0.0 ) &
+          &      antarctica_icesheet_mass_rate_of_change = ( antarctica_icesheet_mass - zantarctica_icesheet_mass_b ) / antarctica_icesheet_timelapsed 
+            antarctica_icesheet_timelapsed = 0.0_wp       
+         ENDIF
+         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass (kg) read in is ', zantarctica_icesheet_mass_in
+         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass (kg) used is    ', antarctica_icesheet_mass
+         IF(lwp) WRITE(numout,*) 'Antarctica icesheet mass rate of change (kg/s) is ', antarctica_icesheet_mass_rate_of_change
+         IF(lwp) WRITE(numout,*) 'Antarctica icesheet seconds lapsed since last change is ', antarctica_icesheet_timelapsed
+      ENDIF
+
       !
       CALL wrk_dealloc( jpi,jpj, ztx, zty, zmsk, zemp, zqns, zqsr, ztx2, zty2 )
@@ -1738 +1894 @@
       REAL(wp), POINTER, DIMENSION(:,:)   ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztmp3, ztmp4   
+
+#if defined key_medusa
+! RSRH Temporarily definitition of array used in medusa coupling until we actually have the medusa code included in the build
+   REAL(wp)  ::  TRC2D_TEST(jpi,jpj)
+#endif
+
       !!----------------------------------------------------------------------
       !
@@ -1879 +2041 @@
       IF( ssnd(jps_co2)%laction )   CALL cpl_snd( jps_co2, isec, RESHAPE ( oce_co2, (/jpi,jpj,1/) ) , info )
       !
+#endif
+
+
+!! JPALM : 03-feb-2015 coupling MEDUSA-UKCA
+!! add soing exactely the same that CO2 fluxes
+!! of PISCES
+!! add CO2-MEDUSA
+!! add DMS-MEDUSA
+!! May add also a coupling MED-UKCA key
+
+
+#if defined key_medusa
+      ! RSRH. We don't want to use magic numbers in the code (i.e. 98 and 221).
+      ! These need moving to a parameter statement (as part of MEDUSA code) or even specifying in a namelist
+      ! so the following code MUST NOT be viewed as anything more than temporary.
+!      IF( ssnd(jps_bio_co2)%laction )   CALL cpl_prism_snd( jps_bio_co2, isec, trc2d(:,:,98:98), info )
+      TRC2D_TEST(:,:) = 0.0
+      IF( ssnd(jps_bio_co2)%laction )   CALL cpl_prism_snd( jps_bio_co2, isec, trc2d_TEST(:,:), info )
+
+      IF( ssnd(jps_bio_dms)%laction )  THEN
+          ! We need to multiply DMS by a conversion factor to get values in the standard units expected in
+          ! the coupling space.
+!         ztmp1(:,: ) = trc2d(:,:,221) * dms_unit_conv
+         ztmp1(:,: ) = trc2d_TEST(:,:) * dms_unit_conv
+         CALL cpl_prism_snd( jps_bio_dms, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
+      ENDIF
 #endif
       !                                                      ! ------------------------- !
@@ -1981 +2169 @@
                                 &         + 0.25 * ( v_ice(ji-1,jj-1) + v_ice(ji,jj-1)                     &
                                 &                  + v_ice(ji-1,jj  ) + v_ice(ji,jj  )  ) *  fr_i(ji,jj)
-                        END DO
+                       END DO
                      END DO
 #if defined key_cice

branches/UKMO/dev_r5518_med_test/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90

@@ -161 +161 @@
       INTEGER, INTENT( in  ) ::   ksbc                ! surface forcing type
       REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
-      REAL(wp) ::   zcoefu, zcoefv, zcoeff            ! local scalar
+      REAL(wp), DIMENSION(:,:,:), POINTER :: ztfrz3d
       INTEGER  ::   ji, jj, jl, jk                    ! dummy loop indices
       !!---------------------------------------------------------------------
@@ -174 +174 @@
       jj_off = INT ( (jpjglo - ny_global) / 2 )
 
-#if defined key_nemocice_decomp
-      ! Pass initial SST from NEMO to CICE so ice is initialised correctly if
-      ! there is no restart file.
-      ! Values from a CICE restart file would overwrite this
-      IF ( .NOT. ln_rstart ) THEN    
-         CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.) 
-      ENDIF  
-#endif
-
-! Initialize CICE
+      ! Initialize CICE
       CALL CICE_Initialize
 
-! Do some CICE consistency checks
+      ! Do some CICE consistency checks
       IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN
          IF ( calc_strair .OR. calc_Tsfc ) THEN
@@ -198 +189 @@
 
 
-! allocate sbc_ice and sbc_cice arrays
-      IF( sbc_ice_alloc()      /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate arrays' )
+      ! allocate sbc_ice and sbc_cice arrays
+      IF( sbc_ice_alloc()      /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_alloc : unable to allocate arrays' )
       IF( sbc_ice_cice_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_cice_alloc : unable to allocate cice arrays' )
 
-! Ensure ocean temperatures are nowhere below freezing if not a NEMO restart
+      ! Ensure that no temperature points are below freezing if not a NEMO restart
       IF( .NOT. ln_rstart ) THEN
-         tsn(:,:,:,jp_tem) = MAX (tsn(:,:,:,jp_tem),Tocnfrz)
+
+         CALL wrk_alloc( jpi,jpj,jpk, ztfrz3d ) 
+         DO jk=1,jpk
+            ztfrz3d(:,:,jk) = eos_fzp( tsn(:,:,jk,jp_sal), fsdept_n(:,:,jk) )
+         ENDDO
+         tsn(:,:,:,jp_tem) = MAX( tsn(:,:,:,jp_tem), ztfrz3d )
          tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
-      ENDIF
-
-      fr_iu(:,:)=0.0
-      fr_iv(:,:)=0.0
+         CALL wrk_dealloc( jpi,jpj,jpk, ztfrz3d ) 
+
+#if defined key_nemocice_decomp
+         ! Pass initial SST from NEMO to CICE so ice is initialised correctly if
+         ! there is no restart file.
+         ! Values from a CICE restart file would overwrite this
+         CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.) 
+#endif
+
+      ENDIF  
+
+      ! calculate surface freezing temperature and send to CICE
+      sstfrz(:,:) = eos_fzp(sss_m(:,:), fsdept_n(:,:,1)) 
+      CALL nemo2cice(sstfrz,Tf, 'T', 1. )
 
       CALL cice2nemo(aice,fr_i, 'T', 1. )
@@ -220 +226 @@
 ! T point to U point
 ! T point to V point
+      fr_iu(:,:)=0.0
+      fr_iv(:,:)=0.0
       DO jj=1,jpjm1
          DO ji=1,jpim1
@@ -348 +356 @@
                ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl)
             ENDDO
-         ELSE
-! emp_ice is set in sbc_cpl_ice_flx as sublimation-snow
-            qla_ice(:,:,1)= - ( emp_ice(:,:)+sprecip(:,:) ) * Lsub
-! End of temporary code
-            DO jj=1,jpj
-               DO ji=1,jpi
-                  IF (fr_i(ji,jj).eq.0.0) THEN
-                     DO jl=1,ncat
-                        ztmpn(ji,jj,jl)=0.0
-                     ENDDO
-                     ! This will then be conserved in CICE
-                     ztmpn(ji,jj,1)=qla_ice(ji,jj,1)
-                  ELSE
-                     DO jl=1,ncat
-                        ztmpn(ji,jj,jl)=qla_ice(ji,jj,1)*a_i(ji,jj,jl)/fr_i(ji,jj)
-                     ENDDO
-                  ENDIF
-               ENDDO
+         ELSE IF (ksbc == jp_purecpl) THEN
+            DO jl=1,ncat
+               ztmpn(:,:,jl)=qla_ice(:,:,jl)*a_i(:,:,jl)
             ENDDO
+    ELSE
+           !In coupled mode - qla_ice calculated in sbc_cpl for each category
+           ztmpn(:,:,1:ncat)=qla_ice(:,:,1:ncat)
          ENDIF
          DO jl=1,ncat

branches/UKMO/dev_r5518_med_test/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

@@ -27 +27 @@
    USE fldread         ! read input field at current time step
 
-
-
    IMPLICIT NONE
    PRIVATE
@@ -90 +88 @@
     INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
     INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
+    REAL(wp)                     ::   zgreenland_fwfisf_sum, zantarctica_fwfisf_sum
     REAL(wp)                     ::   rmin
     REAL(wp)                     ::   zhk
@@ -256 +255 @@
             CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
             fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting) 
+
+            IF( lk_oasis) THEN
+            ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+            IF( ln_coupled_iceshelf_fluxes ) THEN
+
+              ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern
+              ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets
+              ! to preserve total freshwater conservation in coupled models without an active ice sheet model.
+
+               zgreenland_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+               IF( lk_mpp ) CALL mpp_sum( zgreenland_fwfisf_sum )
+               ! use ABS function because we need to preserve the sign of fwfisf
+               WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  &
+              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) &
+              &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) )
+
+               ! check
+               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
+               zgreenland_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+               IF( lk_mpp ) CALL mpp_sum( zgreenland_fwfisf_sum )
+               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum
+
+               zantarctica_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+               IF( lk_mpp ) CALL mpp_sum( zantarctica_fwfisf_sum )
+               ! use ABS function because we need to preserve the sign of fwfisf
+               WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
+              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) &
+              &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) )
+      
+               ! check
+               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
+               zantarctica_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+               IF( lk_mpp ) CALL mpp_sum( zantarctica_fwfisf_sum )
+               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
+
+            ENDIF
+            ENDIF
+
             qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
             stbl(:,:)   = soce
@@ -264 +301 @@
             !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
             fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf
+
+            IF( lk_oasis) THEN
+            ! ln_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
+            IF( ln_coupled_iceshelf_fluxes ) THEN
+
+              ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern
+              ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets
+              ! to preserve total freshwater conservation in coupled models without an active ice sheet model.
+
+               zgreenland_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+               IF( lk_mpp ) CALL mpp_sum( zgreenland_fwfisf_sum )
+               ! use ABS function because we need to preserve the sign of fwfisf
+               WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  &
+              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) &
+              &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) )
+
+               ! check
+               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
+               zgreenland_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
+               IF( lk_mpp ) CALL mpp_sum( zgreenland_fwfisf_sum )
+               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum
+
+               zantarctica_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+               IF( lk_mpp ) CALL mpp_sum( zantarctica_fwfisf_sum )
+               ! use ABS function because we need to preserve the sign of fwfisf
+               WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
+              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) &
+              &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) )
+      
+               ! check
+               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
+               zantarctica_fwfisf_sum = SUM( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
+               IF( lk_mpp ) CALL mpp_sum( zantarctica_fwfisf_sum )
+               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
+
+            ENDIF
+            ENDIF
+
             qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
             !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)              ! heat flux