Changeset 5965 for branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

Timestamp:

2015-12-01T16:35:30+01:00 (8 years ago)

Author:

timgraham

Message:

Upgraded branch to r5518 of trunk (v3.6 stable revision)

File:

• branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90 (modified) (11 diffs)

branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

@@ -9 +9 @@
    !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
    !!             -   !  2010-09  (C. Ethe, G. Madec) Merge TRA-TRC
+   !!            3.6  !  2014-11  (P. Mathiot) isf melting forcing 
    !!----------------------------------------------------------------------
 
@@ -18 +19 @@
    USE dom_oce         ! ocean space domain variables
    USE phycst          ! physical constant
+   USE sbcmod          ! ln_rnf  
+   USE sbcrnf          ! River runoff  
+   USE sbcisf          ! Ice shelf   
    USE traqsr          ! solar radiation penetration
-   USE trdmod_oce      ! ocean trends 
-   USE trdtra          ! ocean trends
+   USE trd_oce         ! trends: ocean variables
+   USE trdtra          ! trends manager: tracers 
+   !
    USE in_out_manager  ! I/O manager
    USE prtctl          ! Print control
-   USE sbcrnf          ! River runoff  
-   USE sbcmod          ! ln_rnf  
    USE iom
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE eosbn2
 
    IMPLICIT NONE
@@ -39 +43 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -91 +95 @@
       !!         where emp, the surface freshwater budget (evaporation minus
       !!         precipitation minus runoff) given in kg/m2/s is divided
-      !!         by rau0 = 1020 kg/m3 (density of sea water) to obtain m/s.    
+      !!         by rau0 (density of sea water) to obtain m/s.    
       !!         Note: even though Fwe does not appear explicitly for 
       !!         temperature in this routine, the heat carried by the water
@@ -107 +111 @@
       !! ** Action  : - Update the 1st level of (ta,sa) with the trend associated
       !!                with the tracer surface boundary condition 
-      !!              - save the trend it in ttrd ('key_trdtra')
+      !!              - send trends to trdtra module (l_trdtra=T)
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
       !!
       INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices  
+      INTEGER  ::   ikt, ikb 
+      INTEGER  ::   nk_isf
       REAL(wp) ::   zfact, z1_e3t, zdep
+      REAL(wp) ::   zalpha, zhk
+      REAL(wp) ::  zt_frz, zpress
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
       !!----------------------------------------------------------------------
@@ -124 +132 @@
       ENDIF
 
-      IF( l_trdtra )   THEN                    !* Save ta and sa trends
+      IF( l_trdtra ) THEN                    !* Save ta and sa trends
          CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
@@ -137 +145 @@
 
       !----------------------------------------
-      !        EMP, EMPS and QNS effects
+      !        EMP, SFX and QNS effects
       !----------------------------------------
       !                                          Set before sbc tracer content fields
@@ -146 +154 @@
               & iom_varid( numror, 'sbc_hc_b', ldstop = .FALSE. ) > 0 ) THEN
             IF(lwp) WRITE(numout,*) '          nit000-1 surface tracer content forcing fields red in the restart file'
-            zfact = 0.5e0
+            zfact = 0.5_wp
             CALL iom_get( numror, jpdom_autoglo, 'sbc_hc_b', sbc_tsc_b(:,:,jp_tem) )   ! before heat content sbc trend
             CALL iom_get( numror, jpdom_autoglo, 'sbc_sc_b', sbc_tsc_b(:,:,jp_sal) )   ! before salt content sbc trend
          ELSE                                         ! No restart or restart not found: Euler forward time stepping
-            zfact = 1.e0
-            sbc_tsc_b(:,:,:) = 0.e0
+            zfact = 1._wp
+            sbc_tsc_b(:,:,:) = 0._wp
          ENDIF
       ELSE                                         ! Swap of forcing fields
          !                                         ! ----------------------
-         zfact = 0.5e0
+         zfact = 0.5_wp
          sbc_tsc_b(:,:,:) = sbc_tsc(:,:,:)
       ENDIF
@@ -182 +190 @@
             END DO
          END DO
-         CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) )                          ! c/d term on sst
-         CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) )                          ! c/d term on sss
+         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * tsn(:,:,1,jp_tem) )   ! c/d term on sst
+         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * tsn(:,:,1,jp_sal) )   ! c/d term on sss
       ENDIF
       ! Concentration dilution effect on (t,s) due to evapouration, precipitation and qns, but not river runoff  
@@ -205 +213 @@
       ENDIF
       !
+      !
+      !----------------------------------------
+      !       Ice Shelf effects (ISF)
+      !     tbl treated as in Losh (2008) JGR
+      !----------------------------------------
+      !
+      IF( nn_isf > 0 ) THEN
+         zfact = 0.5e0
+         DO jj = 2, jpj
+            DO ji = fs_2, fs_jpim1
+         
+               ikt = misfkt(ji,jj)
+               ikb = misfkb(ji,jj)
+   
+               ! level fully include in the ice shelf boundary layer
+               ! if isfdiv, we have to remove heat flux due to inflow at 0oC (as in rnf when you add rnf at sst)
+               ! sign - because fwf sign of evapo (rnf sign of precip)
+               DO jk = ikt, ikb - 1
+               ! compute tfreez for the temperature correction (we add water at freezing temperature)
+!                  zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04
+                  zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress )
+               ! compute trend
+                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                          &
+                     &           + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)          &
+                     &               - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) &
+                     &           * r1_hisf_tbl(ji,jj)
+                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                          &
+                     &           + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj)
+               END DO
+   
+               ! level partially include in ice shelf boundary layer 
+               ! compute tfreez for the temperature correction (we add water at freezing temperature)
+!               zpress = grav*rau0*fsdept(ji,jj,ikb)*1.e-04
+               zt_frz = -1.9 !eos_fzp( tsn(ji,jj,ikb,jp_sal), zpress )
+               ! compute trend
+               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                           &
+                  &              + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)          &
+                  &                  - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) & 
+                  &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)
+               tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal)                                           &
+                  &              + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
+            END DO
+         END DO
+         IF( lrst_oce ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ',   &
+               &                    'at it= ', kt,' date= ', ndastp
+            IF(lwp) WRITE(numout,*) '~~~~'
+            CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:)          )
+            CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem) )
+            CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal) )
+         ENDIF
+      END IF
+      !
       !----------------------------------------
       !        River Runoff effects
@@ -226 +288 @@
       ENDIF
  
-      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
+      IF( l_trdtra )   THEN                      ! send trends for further diagnostics
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
          ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
-         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_nsr, ztrdt )
-         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_nsr, ztrds )
+         CALL trd_tra( kt, 'TRA', jp_tem, jptra_nsr, ztrdt )
+         CALL trd_tra( kt, 'TRA', jp_sal, jptra_nsr, ztrds )
          CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
       ENDIF