Changeset 5708 for branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

Timestamp:

2015-08-26T14:49:08+02:00 (9 years ago)

Author:

davestorkey

Message:

Commit changes for icesheet freshwater input code for coupled models without an active icesheet model.

File:

• branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90 (modified) (5 diffs)

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

@@ -227 +227 @@
          &                  sn_rcv_qns , sn_rcv_emp   , sn_rcv_rnf   , sn_rcv_cal   , sn_rcv_iceflx,  &
          &                  sn_rcv_co2 , sn_rcv_grnm  , sn_rcv_antm  , sn_rcv_ts_ice, nn_cplmodel  ,  &
-         &                  ln_usecplmask
+         &                  ln_usecplmask, ln_coupled_iceshelf_fluxes, rn_greenland_calving_fraction, &
+         &                  rn_antarctica_calving_fraction, rn_iceshelf_fluxes_tolerance
       !!---------------------------------------------------------------------
       !
@@ -285 +286 @@
          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
 
@@ -892 +897 @@
       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 )
       !
@@ -953 +981 @@
       INTEGER  ::   isec                   ! number of seconds since nit000 (assuming rdttra did not change since nit000)
       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
@@ -1237 +1268 @@
       ENDIF
       
-      !                                                        ! land ice masses
-      IF( srcv(jpr_grnm)%laction )   greenland_mass_array(:,:) = frcv(jpr_grnm)%z3(:,:,1)
-      IF( srcv(jpr_antm)%laction )   antarctica_mass_array(:,:) = frcv(jpr_antm)%z3(:,:,1)
-      greenland_mass = greenland_mass_array(1,1)
-      antarctica_mass = antarctica_mass_array(1,1)
-         
-      WRITE(numout,*) 'Greenland mass is ', greenland_mass
-      WRITE(numout,*) 'Antarctica mass is ', antarctica_mass
+      !                                                        ! 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 )/( greenland_icesheet_mass + zepsilon ) > zepsilon ) THEN
+            zgreenland_icesheet_mass_b = greenland_icesheet_mass
+            greenland_icesheet_mass = zgreenland_icesheet_mass_in
+            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) 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 )/( antarctica_icesheet_mass + zepsilon ) > zepsilon ) THEN
+            zantarctica_icesheet_mass_b = antarctica_icesheet_mass
+            antarctica_icesheet_mass = zantarctica_icesheet_mass_in
+            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) 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 )