Changeset 5663 for branches/UKMO/dev_r5518_CICE_coupling_GSI7_GSI8/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90

Timestamp:

2015-08-03T15:53:06+02:00 (9 years ago)

Author:

dancopsey

Message:

Merged in Alex West's GSI8 changes from eld259:/data/local/hadax/FCM_working/NEMO/Multilayers/NEMO3.6_stable/UKMO1_CICE_coupling_GSI7_GSI8

File:

• branches/UKMO/dev_r5518_CICE_coupling_GSI7_GSI8/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90 (modified) (13 diffs)

branches/UKMO/dev_r5518_CICE_coupling_GSI7_GSI8/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,                    &
                 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 +171 @@
       INTEGER, INTENT( in  ) ::   ksbc                ! surface forcing type
       REAL(wp), DIMENSION(:,:), POINTER :: ztmp1, ztmp2
+      REAL(wp), DIMENSION(:,:,:), POINTER :: ztfrz3d
       REAL(wp) ::   zcoefu, zcoefv, zcoeff            ! local scalar
       INTEGER  ::   ji, jj, jl, jk                    ! dummy loop indices
@@ -173 +184 @@
       ji_off = INT ( (jpiglo - nx_global) / 2 )
       jj_off = INT ( (jpjglo - ny_global) / 2 )
+
+      ! Set freezing temperatures and ensure consistencey between NEMO and CICE
+      CALL wrk_alloc( jpi,jpj,jpk, ztfrz3d ) 
+      DO jk=1,jpk
+         ztfrz3d(:,:,jk) = eos_fzp( tsn(:,:,jk,jp_sal), fsdept_n(:,:,jk) )
+      ENDDO
+
+      !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),ztfrz3d)
+         tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
+      ENDIF
 
 #if defined key_nemocice_decomp
@@ -202 +225 @@
       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
-      IF( .NOT. ln_rstart ) THEN
-         tsn(:,:,:,jp_tem) = MAX (tsn(:,:,:,jp_tem),Tocnfrz)
-         tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)
-      ENDIF
+      ! Populate the surface freezing temperature array
+      sstfrz(:,:)=ztfrz3d(:,:,1)
 
       fr_iu(:,:)=0.0
@@ -283 +303 @@
  
       CALL wrk_dealloc( jpi,jpj, ztmp1, ztmp2 )
+      CALL wrk_dealloc( jpi,jpj,jpk, ztfrz3d ) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('cice_sbc_init')
@@ -343 +364 @@
          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
+         IF (ksbc == jp_flx .OR. ksbc == jp_purecpl) THEN
             DO jl=1,ncat
                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
-            ENDDO
+           !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 +391 @@
 ! 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 +400 @@
 ! 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
@@ -442 +460 @@
       CALL nemo2cice(ztmp,frain,'T', 1. ) 
 
+! Recalculate freezing temperature and send to CICE 
+      sstfrz(:,:)=eos_fzp(sss_m(:,:), 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 +473 @@
       CALL nemo2cice(sst_m,sst,'T', 1. )
       CALL nemo2cice(sss_m,sss,'T', 1. )
+
+! Sea ice surface skin temperature
+      DO jl=1,ncat
+        CALL nemo2cice(tsfc_ice(:,:,jl), trcrn(:,:,nt_tsfc,jl,:),'T',1.)
+      ENDDO 
 
 ! x comp and y comp of surface ocean current
@@ -730 +755 @@
          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')