Changeset 9987 for branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90

Timestamp:

2018-07-23T11:33:03+02:00 (6 years ago)

Author:

emmafiedler

Message:

Merge with GO6 FOAMv14 package branch r9288

File:

• branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 (modified) (18 diffs)

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

@@ -15 +15 @@
    USE dom_oce         ! ocean space and time domain
    USE domvvl
-   USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic
+   USE eosbn2, only : eos_fzp ! Function to calculate freezing point of seawater
+   USE phycst, only : rcp, rau0, r1_rau0, rhosn, rhoic, rt0
    USE in_out_manager  ! I/O manager
    USE iom, ONLY : iom_put,iom_use              ! I/O manager library !!Joakim edit
@@ -37 +38 @@
    USE ice_gather_scatter
    USE ice_calendar, only: dt
+# if defined key_cice4
    USE ice_state, only: aice,aicen,uvel,vvel,vsno,vsnon,vice,vicen
-# if defined key_cice4
    USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow,  &
                 strocnxT,strocnyT,                               & 
@@ -45 +46 @@
                 flatn_f,fsurfn_f,fcondtopn_f,                    &
                 uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl,   &
-                swvdr,swvdf,swidr,swidf
+                swvdr,swvdf,swidr,swidf,Tf
    USE ice_therm_vertical, only: calc_Tsfc
 #else
+   USE ice_state, only: aice,aicen,uvel,nt_hpnd,trcrn,vvel,vsno,&
+                vsnon,vice,vicen,nt_Tsfc
    USE ice_flux, only: strax,stray,strocnx,strocny,frain,fsnow,  &
                 strocnxT,strocnyT,                               & 
-                sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai,     &
-                fresh_ai,fhocn_ai,fswthru_ai,frzmlt,          &
+                sst,sss,uocn,vocn,ss_tltx,ss_tlty,fsalt_ai,      &
+                fresh_ai,fhocn_ai,fswthru_ai,frzmlt,             &
                 flatn_f,fsurfn_f,fcondtopn_f,                    &
+#ifdef key_asminc
+                daice_da,fresh_da,fsalt_da,                    &
+#endif
                 uatm,vatm,wind,fsw,flw,Tair,potT,Qa,rhoa,zlvl,   &
-                swvdr,swvdf,swidr,swidf
-   USE ice_therm_shared, only: calc_Tsfc
+                swvdr,swvdf,swidr,swidf,Tf,                      &
+      !! When using NEMO with CICE, this change requires use of 
+      !! one of the following two CICE branches:
+      !! - at CICE5.0,   hadax/r1015_GSI8_with_GSI7
+      !! - at CICE5.1.2, hadax/vn5.1.2_GSI8
+                keffn_top,Tn_top
+
+   USE ice_therm_shared, only: calc_Tsfc, heat_capacity
+   USE ice_shortwave, only: apeffn
 #endif
    USE ice_forcing, only: frcvdr,frcvdf,frcidr,frcidf
@@ -161 +174 @@
       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 +187 @@
       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 +202 @@
 
 
-! 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
+             CALL eos_fzp( tsn(:,:,jk,jp_sal), ztfrz3d(:,:,jk), 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
+      CALL  eos_fzp(sss_m(:,:), sstfrz(:,:), fsdept_n(:,:,1)) 
+      CALL nemo2cice(sstfrz,Tf, 'T', 1. )
 
       CALL cice2nemo(aice,fr_i, 'T', 1. )
@@ -220 +239 @@
 ! 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
@@ -283 +304 @@
  
       CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
-      !
+
+#if defined key_asminc
+      ! Initialize fresh water and salt fluxes from data assim   
+      !  and data assimilation index to cice 
+      nfresh_da(:,:) = 0.0   
+      nfsalt_da(:,:) = 0.0   
+      ndaice_da(:,:) = 0.0         
+#endif
+      !
+      ! In coupled mode get extra fields from CICE for passing back to atmosphere
+ 
+      IF ( ksbc == jp_purecpl ) CALL cice_sbc_hadgam(nit000)
+      ! 
       IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_init')
       !
@@ -343 +376 @@
          CALL nemo2cice(ztmp,stray,'F', -1. )
 
+
+! Alex West: From configuration GSI8 onwards, when NEMO is used with CICE in
+! HadGEM3 the 'time-travelling ice' coupling approach is used, whereby 
+! atmosphere-ice fluxes are passed as pseudo-local values, formed by dividing
+! gridbox mean fluxes in the UM by future ice concentration obtained through  
+! OASIS.  This allows for a much more realistic apportionment of energy through
+! the ice - and conserves energy.
+! Therefore the fluxes are now divided by ice concentration in the coupled
+! formulation (jp_purecpl) as well as for jp_flx.  This NEMO branch should only
+! be used at UM10.2 onwards (unless an explicit GSI8 UM branch is included), at
+! which point the GSI8 UM changes were committed.
+
 ! Surface downward latent heat flux (CI_5)
          IF (ksbc == jp_flx) THEN
@@ -348 +393 @@
                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
             CALL nemo2cice(ztmpn(:,:,jl),flatn_f(:,:,jl,:),'T', 1. )
@@ -373 +407 @@
 ! GBM conductive flux through ice (CI_6)
 !  Convert to GBM
-            IF (ksbc == jp_flx) THEN
+            IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
                ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl)
             ELSE
@@ -382 +416 @@
 ! GBM surface heat flux (CI_7)
 !  Convert to GBM
-            IF (ksbc == jp_flx) THEN
+            IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
                ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl) 
             ELSE
@@ -431 +465 @@
       ENDIF
 
+#if defined key_asminc
+!Ice concentration change (from assimilation)
+      ztmp(:,:)=ndaice_da(:,:)*tmask(:,:,1)
+      Call nemo2cice(ztmp,daice_da,'T', 1. )
+#endif 
+
 ! Snowfall
 ! Ensure fsnow is positive (as in CICE routine prepare_forcing)
       IF( iom_use('snowpre') )   CALL iom_put('snowpre',MAX( (1.0-fr_i(:,:))*sprecip(:,:) ,0.0)) !!Joakim edit  
-      ztmp(:,:)=MAX(fr_i(:,:)*sprecip(:,:),0.0)  
+      IF( kt == nit000 .AND. lwp )  THEN
+         WRITE(numout,*) 'sprecip weight, rn_sfac=', rn_sfac
+      ENDIF
+      ztmp(:,:)=MAX(fr_i(:,:)*rn_sfac*sprecip(:,:),0.0)  
       CALL nemo2cice(ztmp,fsnow,'T', 1. ) 
 
@@ -442 +485 @@
       CALL nemo2cice(ztmp,frain,'T', 1. ) 
 
+! Recalculate freezing temperature and send to CICE 
+      CALL eos_fzp(sss_m(:,:), sstfrz(:,:), fsdept_n(:,:,1)) 
+      CALL nemo2cice(sstfrz,Tf,'T', 1. )
+
 ! Freezing/melting potential
 ! Calculated over NEMO leapfrog timestep (hence 2*dt)
-      nfrzmlt(:,:)=rau0*rcp*fse3t_m(:,:)*(Tocnfrz-sst_m(:,:))/(2.0*dt)
-
-      ztmp(:,:) = nfrzmlt(:,:)
-      CALL nemo2cice(ztmp,frzmlt,'T', 1. )
+      nfrzmlt(:,:)=rau0*rcp*fse3t_m(:,:)*(sstfrz(:,:)-sst_m(:,:))/(2.0*dt) 
+      CALL nemo2cice(nfrzmlt,frzmlt,'T', 1. )
 
 ! SST  and SSS
@@ -453 +498 @@
       CALL nemo2cice(sst_m,sst,'T', 1. )
       CALL nemo2cice(sss_m,sss,'T', 1. )
+
+      IF( ksbc == jp_purecpl ) THEN
+! Sea ice surface skin temperature
+         DO jl=1,ncat
+           CALL nemo2cice(tsfc_ice(:,:,jl), trcrn(:,:,nt_tsfc,jl,:),'T',1.)
+         ENDDO 
+      ENDIF
 
 ! x comp and y comp of surface ocean current
@@ -685 +737 @@
       ENDIF
 
+#if defined key_asminc
+! Import fresh water and salt flux due to seaice da
+      CALL cice2nemo(fresh_da, nfresh_da,'T',1.0)
+      CALL cice2nemo(fsalt_da, nfsalt_da,'T',1.0)
+#endif
+
 ! Release work space
 
@@ -708 +766 @@
       IF( nn_timing == 1 )  CALL timing_start('cice_sbc_hadgam')
       !
-      IF( kt == nit000 )  THEN
-         IF(lwp) WRITE(numout,*)'cice_sbc_hadgam'
-         IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
-      ENDIF
-
       !                                         ! =========================== !
       !                                         !   Prepare Coupling fields   !
@@ -730 +783 @@
          CALL cice2nemo(vicen(:,:,jl,:),ht_i(:,:,jl),'T', 1. )
       ENDDO
+
+#if ! defined key_cice4
+! Meltpond fraction and depth
+      DO jl = 1,ncat
+         CALL cice2nemo(apeffn(:,:,jl,:),a_p(:,:,jl),'T', 1. )
+         CALL cice2nemo(trcrn(:,:,nt_hpnd,jl,:),ht_p(:,:,jl),'T', 1. )
+      ENDDO
+#endif
+
+
+! If using multilayers thermodynamics in CICE then get top layer temperature
+! and effective conductivity       
+!! When using NEMO with CICE, this change requires use of 
+!! one of the following two CICE branches:
+!! - at CICE5.0,   hadax/r1015_GSI8_with_GSI7
+!! - at CICE5.1.2, hadax/vn5.1.2_GSI8
+      IF (heat_capacity) THEN
+         DO jl = 1,ncat
+            CALL cice2nemo(Tn_top(:,:,jl,:),tn_ice(:,:,jl),'T', 1. )
+            CALL cice2nemo(keffn_top(:,:,jl,:),kn_ice(:,:,jl),'T', 1. )
+         ENDDO
+! Convert surface temperature to Kelvin
+         tn_ice(:,:,:)=tn_ice(:,:,:)+rt0
+      ELSE
+         tn_ice(:,:,:) = 0.0
+         kn_ice(:,:,:) = 0.0
+      ENDIF       
+
       !
       IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_hadgam')