New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 5682 for branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Ignore:
Timestamp:
2015-08-12T17:46:45+02:00 (9 years ago)
Author:
mattmartin
Message:

OBS simplification changes committed to branch after running SETTE tests to make sure we get the same results as the trunk for ORCA2_LIM_OBS.

Location:
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA
Files:
25 edited

Legend:

Unmodified
Added
Removed

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r4990 r5682  
    2222   !!             -   ! 2013-04  (F. Roquet, G. Madec)  add eos_rab, change bn2 computation and reorganize the module 
    2323   !!             -   ! 2014-09  (F. Roquet)  add TEOS-10, S-EOS, and modify EOS-80 
     24   !!             -   ! 2015-06  (P.A. Bouttier) eos_fzp functions changed to subroutines for AGRIF 
    2425   !!---------------------------------------------------------------------- 
    2526 
    … …  
    4748   USE lbclnk         ! ocean lateral boundary conditions 
    4849   USE timing          ! Timing 
     50   USE stopar          ! Stochastic T/S fluctuations 
     51   USE stopts          ! Stochastic T/S fluctuations 
    4952 
    5053   IMPLICIT NONE 
    … …  
    7275   PUBLIC   eos_init       ! called by istate module 
    7376 
    74    !                                          !!* Namelist (nameos) * 
    75    INTEGER , PUBLIC ::   nn_eos   = 0         !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 
    76    LOGICAL , PUBLIC ::   ln_useCT  = .FALSE. ! determine if eos_pt_from_ct is used to compute sst_m 
     77   !                                !!* Namelist (nameos) * 
     78   INTEGER , PUBLIC ::   nn_eos     ! = 0/1/2 type of eq. of state and Brunt-Vaisala frequ. 
     79   LOGICAL , PUBLIC ::   ln_useCT  ! determine if eos_pt_from_ct is used to compute sst_m 
    7780 
    7881   !                                   !!!  simplified eos coefficients 
    … …  
    313316      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pdep   ! depth                      [m] 
    314317      ! 
    315       INTEGER  ::   ji, jj, jk                ! dummy loop indices 
    316       REAL(wp) ::   zt , zh , zs , ztm        ! local scalars 
    317       REAL(wp) ::   zn , zn0, zn1, zn2, zn3   !   -      - 
     318      INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices 
     319      INTEGER  ::   jdof 
     320      REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars 
     321      REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      - 
     322      REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors 
    318323      !!---------------------------------------------------------------------- 
    319324      ! 
    … …  
    324329      CASE( -1, 0 )                !==  polynomial TEOS-10 / EOS-80 ==! 
    325330         ! 
    326          DO jk = 1, jpkm1 
    327             DO jj = 1, jpj 
    328                DO ji = 1, jpi 
    329                   ! 
    330                   zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
    331                   zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
    332                   zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
    333                   ztm = tmask(ji,jj,jk)                                         ! tmask 
    334                   ! 
    335                   zn3 = EOS013*zt   & 
    336                      &   + EOS103*zs+EOS003 
    337                      ! 
    338                   zn2 = (EOS022*zt   & 
    339                      &   + EOS112*zs+EOS012)*zt   & 
    340                      &   + (EOS202*zs+EOS102)*zs+EOS002 
    341                      ! 
    342                   zn1 = (((EOS041*zt   & 
    343                      &   + EOS131*zs+EOS031)*zt   & 
    344                      &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
    345                      &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
    346                      &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
    347                      ! 
    348                   zn0 = (((((EOS060*zt   & 
    349                      &   + EOS150*zs+EOS050)*zt   & 
    350                      &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
    351                      &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
    352                      &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
    353                      &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
    354                      &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
    355                      ! 
    356                   zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
    357                   ! 
    358                   prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface 
    359                   ! 
    360                   prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm      ! density anomaly (masked) 
     331         ! Stochastic equation of state 
     332         IF ( ln_sto_eos ) THEN 
     333            ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 
     334            ALLOCATE(zn_sto(1:2*nn_sto_eos)) 
     335            ALLOCATE(zsign(1:2*nn_sto_eos)) 
     336            DO jsmp = 1, 2*nn_sto_eos, 2 
     337              zsign(jsmp)   = 1._wp 
     338              zsign(jsmp+1) = -1._wp 
     339            END DO 
     340            ! 
     341            DO jk = 1, jpkm1 
     342               DO jj = 1, jpj 
     343                  DO ji = 1, jpi 
     344                     ! 
     345                     ! compute density (2*nn_sto_eos) times: 
     346                     ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 
     347                     ! (2) for t-dt, s-ds (with the opposite fluctuation) 
     348                     DO jsmp = 1, nn_sto_eos*2 
     349                        jdof   = (jsmp + 1) / 2 
     350                        zh     = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     351                        zt     = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0    ! temperature 
     352                        zstemp = pts  (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 
     353                        zs     = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 )   ! square root salinity 
     354                        ztm    = tmask(ji,jj,jk)                                         ! tmask 
     355                        ! 
     356                        zn3 = EOS013*zt   & 
     357                           &   + EOS103*zs+EOS003 
     358                           ! 
     359                        zn2 = (EOS022*zt   & 
     360                           &   + EOS112*zs+EOS012)*zt   & 
     361                           &   + (EOS202*zs+EOS102)*zs+EOS002 
     362                           ! 
     363                        zn1 = (((EOS041*zt   & 
     364                           &   + EOS131*zs+EOS031)*zt   & 
     365                           &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     366                           &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     367                           &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     368                           ! 
     369                        zn0_sto(jsmp) = (((((EOS060*zt   & 
     370                           &   + EOS150*zs+EOS050)*zt   & 
     371                           &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     372                           &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     373                           &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     374                           &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     375                           &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     376                           ! 
     377                        zn_sto(jsmp)  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp) 
     378                     END DO 
     379                     ! 
     380                     ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 
     381                     prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 
     382                     DO jsmp = 1, nn_sto_eos*2 
     383                        prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp)                      ! potential density referenced at the surface 
     384                        ! 
     385                        prd(ji,jj,jk) = prd(ji,jj,jk) + (  zn_sto(jsmp) * r1_rau0 - 1._wp  )   ! density anomaly (masked) 
     386                     END DO 
     387                     prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 
     388                     prd  (ji,jj,jk) = 0.5_wp * prd  (ji,jj,jk) * ztm / nn_sto_eos 
     389                  END DO 
    361390               END DO 
    362391            END DO 
    363          END DO 
    364          ! 
     392            DEALLOCATE(zn0_sto,zn_sto,zsign) 
     393         ! Non-stochastic equation of state 
     394         ELSE 
     395            DO jk = 1, jpkm1 
     396               DO jj = 1, jpj 
     397                  DO ji = 1, jpi 
     398                     ! 
     399                     zh  = pdep(ji,jj,jk) * r1_Z0                                  ! depth 
     400                     zt  = pts (ji,jj,jk,jp_tem) * r1_T0                           ! temperature 
     401                     zs  = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity 
     402                     ztm = tmask(ji,jj,jk)                                         ! tmask 
     403                     ! 
     404                     zn3 = EOS013*zt   & 
     405                        &   + EOS103*zs+EOS003 
     406                        ! 
     407                     zn2 = (EOS022*zt   & 
     408                        &   + EOS112*zs+EOS012)*zt   & 
     409                        &   + (EOS202*zs+EOS102)*zs+EOS002 
     410                        ! 
     411                     zn1 = (((EOS041*zt   & 
     412                        &   + EOS131*zs+EOS031)*zt   & 
     413                        &   + (EOS221*zs+EOS121)*zs+EOS021)*zt   & 
     414                        &   + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt   & 
     415                        &   + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 
     416                        ! 
     417                     zn0 = (((((EOS060*zt   & 
     418                        &   + EOS150*zs+EOS050)*zt   & 
     419                        &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   & 
     420                        &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   & 
     421                        &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   & 
     422                        &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   & 
     423                        &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 
     424                        ! 
     425                     zn  = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 
     426                     ! 
     427                     prhop(ji,jj,jk) = zn0 * ztm                           ! potential density referenced at the surface 
     428                     ! 
     429                     prd(ji,jj,jk) = (  zn * r1_rau0 - 1._wp  ) * ztm      ! density anomaly (masked) 
     430                  END DO 
     431               END DO 
     432            END DO 
     433         ENDIF 
     434          
    365435      CASE( 1 )                !==  simplified EOS  ==! 
    366436         ! 
    … …  
    922992 
    923993 
    924    FUNCTION eos_fzp_2d( psal, pdep ) RESULT( ptf ) 
     994   SUBROUTINE  eos_fzp_2d( psal, ptf, pdep ) 
    925995      !!---------------------------------------------------------------------- 
    926996      !!                 ***  ROUTINE eos_fzp  *** 
    … …  
    9361006      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )           ::   psal   ! salinity   [psu] 
    9371007      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ), OPTIONAL ::   pdep   ! depth      [m] 
    938       REAL(wp), DIMENSION(jpi,jpj)                          ::   ptf   ! freezing temperature [Celcius] 
     1008      REAL(wp), DIMENSION(jpi,jpj), INTENT(out  )           ::   ptf    ! freezing temperature [Celcius] 
    9391009      ! 
    9401010      INTEGER  ::   ji, jj   ! dummy loop indices 
    … …  
    9691039         nstop = nstop + 1 
    9701040         ! 
    971       END SELECT 
    972       ! 
    973    END FUNCTION eos_fzp_2d 
    974  
    975   FUNCTION eos_fzp_0d( psal, pdep ) RESULT( ptf ) 
     1041      END SELECT       
     1042      ! 
     1043  END SUBROUTINE eos_fzp_2d 
     1044 
     1045  SUBROUTINE eos_fzp_0d( psal, ptf, pdep ) 
    9761046      !!---------------------------------------------------------------------- 
    9771047      !!                 ***  ROUTINE eos_fzp  *** 
    … …  
    9851055      !! Reference  :   UNESCO tech. papers in the marine science no. 28. 1978 
    9861056      !!---------------------------------------------------------------------- 
    987       REAL(wp), INTENT(in)           ::   psal   ! salinity   [psu] 
    988       REAL(wp), INTENT(in), OPTIONAL ::   pdep   ! depth      [m] 
    989       REAL(wp)                       ::   ptf   ! freezing temperature [Celcius] 
     1057      REAL(wp), INTENT(in )           ::   psal         ! salinity   [psu] 
     1058      REAL(wp), INTENT(in ), OPTIONAL ::   pdep         ! depth      [m] 
     1059      REAL(wp), INTENT(out)           ::   ptf          ! freezing temperature [Celcius] 
    9901060      ! 
    9911061      REAL(wp) :: zs   ! local scalars 
    … …  
    10171087      END SELECT 
    10181088      ! 
    1019    END FUNCTION eos_fzp_0d 
     1089   END SUBROUTINE eos_fzp_0d 
    10201090 
    10211091 
    … …  
    11831253            WRITE(numout,*) '             model uses Conservative Temperature' 
    11841254            WRITE(numout,*) '             Important: model must be initialized with CT and SA fields' 
     1255         ELSE 
     1256            WRITE(numout,*) '             model does not use Conservative Temperature' 
    11851257         ENDIF 
    11861258      ENDIF 
    … …  
    15891661      END SELECT 
    15901662      ! 
     1663      rau0_rcp    = rau0 * rcp  
    15911664      r1_rau0     = 1._wp / rau0 
    15921665      r1_rcp      = 1._wp / rcp 
    1593       r1_rau0_rcp = 1._wp / ( rau0 * rcp ) 
     1666      r1_rau0_rcp = 1._wp / rau0_rcp  
    15941667      ! 
    15951668      IF(lwp) WRITE(numout,*) 
    … …  
    15971670      IF(lwp) WRITE(numout,*) '          1. / rau0                        r1_rau0  = ', r1_rau0, ' m^3/kg' 
    15981671      IF(lwp) WRITE(numout,*) '          ocean specific heat                 rcp   = ', rcp    , ' J/Kelvin' 
     1672      IF(lwp) WRITE(numout,*) '          rau0 * rcp                       rau0_rcp = ', rau0_rcp 
    15991673      IF(lwp) WRITE(numout,*) '          1. / ( rau0 * rcp )           r1_rau0_rcp = ', r1_rau0_rcp 
    16001674      ! 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

    r4990 r5682  
    2626   USE cla             ! cross land advection      (cla_traadv     routine) 
    2727   USE ldftra_oce      ! lateral diffusion coefficient on tracers 
     28   ! 
    2829   USE in_out_manager  ! I/O manager 
    2930   USE iom             ! I/O module 
    … …  
    3334   USE timing          ! Timing 
    3435   USE sbc_oce 
     36   USE diaptr          ! Poleward heat transport  
    3537 
    3638 
    … …  
    111113      ! 
    112114      IF( ln_mle    )   CALL tra_adv_mle( kt, nit000, zun, zvn, zwn, 'TRA' )    ! add the mle transport (if necessary) 
     115      ! 
    113116      CALL iom_put( "uocetr_eff", zun )                                         ! output effective transport       
    114117      CALL iom_put( "vocetr_eff", zvn ) 
    115118      CALL iom_put( "wocetr_eff", zwn ) 
    116  
     119      ! 
     120      IF( ln_diaptr )   CALL dia_ptr( zvn )                                     ! diagnose the effective MSF  
     121      ! 
     122    
    117123      SELECT CASE ( nadv )                            !==  compute advection trend and add it to general trend  ==! 
    118       CASE ( 1 )   ;    CALL tra_adv_cen2  ( kt, nit000, 'TRA',         zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  2nd order centered 
    119       CASE ( 2 )   ;    CALL tra_adv_tvd   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD  
    120       CASE ( 3 )   ;    CALL tra_adv_muscl ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb,      tsa, jpts, ln_traadv_msc_ups )   !  MUSCL  
    121       CASE ( 4 )   ;    CALL tra_adv_muscl2( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  MUSCL2  
    122       CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  UBS  
    123       CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  QUICKEST  
    124       CASE ( 7 )   ;   CALL tra_adv_tvd_zts( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD ZTS 
     124      CASE ( 1 )   ;    CALL tra_adv_cen2   ( kt, nit000, 'TRA',         zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  2nd order centered 
     125      CASE ( 2 )   ;    CALL tra_adv_tvd    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD  
     126      CASE ( 3 )   ;    CALL tra_adv_muscl  ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb,      tsa, jpts, ln_traadv_msc_ups )   !  MUSCL  
     127      CASE ( 4 )   ;    CALL tra_adv_muscl2 ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  MUSCL2  
     128      CASE ( 5 )   ;    CALL tra_adv_ubs    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  UBS  
     129      CASE ( 6 )   ;    CALL tra_adv_qck    ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  QUICKEST  
     130      CASE ( 7 )   ;    CALL tra_adv_tvd_zts( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts )   !  TVD ZTS 
    125131      ! 
    126132      CASE (-1 )                                      !==  esopa: test all possibility with control print  ==! 
    … …  
    206212      IF( lk_esopa         )   ioptio =          1 
    207213 
    208       IF( ( ln_traadv_muscl .OR. ln_traadv_muscl2 .OR. ln_traadv_ubs .OR. ln_traadv_qck ) .AND. nn_isf .NE. 0 )  & 
    209       &   CALL ctl_stop( 'Only traadv_cen2 and traadv_tvd is compatible with ice shelf cavity') 
     214      IF( ( ln_traadv_muscl .OR. ln_traadv_muscl2 .OR. ln_traadv_ubs .OR. ln_traadv_qck .OR. ln_traadv_tvd_zts )   & 
     215         .AND. ln_isfcav )   CALL ctl_stop( 'Only traadv_cen2 and traadv_tvd is compatible with ice shelf cavity') 
    210216 
    211217      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namtra_adv' ) 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90

    r4990 r5682  
    173173         END DO  
    174174      END DO  
    175       zfzp(:,:) = eos_fzp( tsn(:,:,1,jp_sal), zpres(:,:) ) 
     175      CALL eos_fzp( tsn(:,:,1,jp_sal), zfzp(:,:), zpres(:,:) ) 
    176176      DO jk = 1, jpk 
    177177         DO jj = 1, jpj 
    … …  
    279279         END IF 
    280280         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    281          IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN   
    282            IF( jn == jp_tem )   htr_adv(:) = ptr_vj( zwy(:,:,:) ) 
    283            IF( jn == jp_sal )   str_adv(:) = ptr_vj( zwy(:,:,:) ) 
     281         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
     282           IF( jn == jp_tem )   htr_adv(:) = ptr_sj( zwy(:,:,:) ) 
     283           IF( jn == jp_sal )   str_adv(:) = ptr_sj( zwy(:,:,:) ) 
    284284         ENDIF 
    285285         ! 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_eiv.F90

    r4990 r5682  
    212212      CHARACTER(len=3) ::   cdtype 
    213213      REAL, DIMENSION(:,:,:) ::   pun, pvn, pwn 
    214       WRITE(*,*) 'tra_adv_eiv: You should not have seen this print! error?', kt, cdtype, pun(1,1,1), pvn(1,1,1), pwn(1,1,1) 
     214      WRITE(*,*) 'tra_adv_eiv: You should not have seen this print! error?',   & 
     215         &    kt, cdtype, pun(1,1,1), pvn(1,1,1), pwn(1,1,1) 
    215216   END SUBROUTINE tra_adv_eiv 
    216217#endif 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mle.F90

    • Property svn:keywords set to Id
    r4835 r5682  
    5353   !!---------------------------------------------------------------------- 
    5454   !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
    55    !! $Id:$ 
     55   !! $Id$ 
    5656   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    5757   !!---------------------------------------------------------------------- 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90

    r4990 r5682  
    2121   USE trdtra         ! tracers trends manager 
    2222   USE dynspg_oce     ! choice/control of key cpp for surface pressure gradient 
    23    USE sbcrnf          ! river runoffs 
     23   USE sbcrnf         ! river runoffs 
    2424   USE diaptr         ! poleward transport diagnostics 
    2525   ! 
    … …  
    219219         END IF 
    220220         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    221          IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN   
    222             IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) ) 
    223             IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) ) 
     221         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
     222            IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) 
     223            IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) 
    224224         ENDIF 
    225225 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90

    r4990 r5682  
    200200 
    201201         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    202          IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
    203             IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) ) 
    204             IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) ) 
     202         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
     203            IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) 
     204            IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) 
    205205         ENDIF 
    206206 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90

    r4990 r5682  
    355355         IF( l_trd )   CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) ) 
    356356         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    357          IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN   
    358            IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) ) 
    359            IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) ) 
     357         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
     358           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) 
     359           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) 
    360360         ENDIF 
    361361         ! 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90

    r4990 r5682  
    106106      ENDIF 
    107107      ! 
    108       zwi(:,:,:) = 0.e0 ; zwz(:,:,:) = 0.e0 
     108      zwi(:,:,:) = 0.e0 ;  
    109109      ! 
    110110      !                                                          ! =========== 
    111111      DO jn = 1, kjpt                                            ! tracer loop 
    112112         !                                                       ! =========== 
    113          ! 1. Bottom value : flux set to zero 
     113         ! 1. Bottom and k=1 value : flux set to zero 
    114114         ! ---------------------------------- 
    115115         zwx(:,:,jpk) = 0.e0    ;    zwz(:,:,jpk) = 0.e0 
    116116         zwy(:,:,jpk) = 0.e0    ;    zwi(:,:,jpk) = 0.e0 
    117  
     117           
     118         zwz(:,:,1  ) = 0._wp 
    118119         ! 2. upstream advection with initial mass fluxes & intermediate update 
    119120         ! -------------------------------------------------------------------- 
    … …  
    134135 
    135136         ! upstream tracer flux in the k direction 
     137         ! Interior value 
     138         DO jk = 2, jpkm1 
     139            DO jj = 1, jpj 
     140               DO ji = 1, jpi 
     141                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
     142                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
     143                  zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 
     144               END DO 
     145            END DO 
     146         END DO 
    136147         ! Surface value 
    137148         IF( lk_vvl ) THEN    
    138             DO jj = 1, jpj 
    139                DO ji = 1, jpi 
    140                   zwz(ji,jj, mikt(ji,jj) ) = 0.e0                         ! volume variable 
    141                END DO 
    142             END DO 
     149            IF ( ln_isfcav ) THEN 
     150               DO jj = 1, jpj 
     151                  DO ji = 1, jpi 
     152                     zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable 
     153                  END DO 
     154               END DO 
     155            ELSE 
     156               zwz(:,:,1) = 0.e0          ! volume variable 
     157            END IF 
    143158         ELSE                 
    144             DO jj = 1, jpj 
    145                DO ji = 1, jpi 
    146                   zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
    147                END DO 
    148             END DO    
     159            IF ( ln_isfcav ) THEN 
     160               DO jj = 1, jpj 
     161                  DO ji = 1, jpi 
     162                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
     163                  END DO 
     164               END DO    
     165            ELSE 
     166               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface 
     167            END IF 
    149168         ENDIF 
    150          ! Interior value 
    151          DO jj = 1, jpj 
    152             DO ji = 1, jpi 
    153                DO jk = mikt(ji,jj)+1, jpkm1 
    154                   zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
    155                   zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
    156                   zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) 
    157                END DO 
    158             END DO 
    159          END DO 
    160169 
    161170         ! total advective trend 
    … …  
    184193         END IF 
    185194         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    186          IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN   
    187            IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) ) 
    188            IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) ) 
     195         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
     196           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) 
     197           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) 
    189198         ENDIF 
    190199 
    … …  
    202211       
    203212         ! antidiffusive flux on k 
    204          zwz(:,:,1) = 0.e0         ! Surface value 
    205          ! 
    206          DO jj = 1, jpj 
    207             DO ji = 1, jpi 
    208                ik=mikt(ji,jj) 
    209                ! surface value 
    210                zwz(ji,jj,1:ik) = 0.e0 
    211                ! Interior value 
    212                DO jk = mikt(ji,jj)+1, jpkm1                     
     213         ! Interior value 
     214         DO jk = 2, jpkm1                     
     215            DO jj = 1, jpj 
     216               DO ji = 1, jpi 
    213217                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk) 
    214218               END DO 
    215219            END DO 
    216220         END DO 
     221         ! surface value 
     222         IF ( ln_isfcav ) THEN 
     223            DO jj = 1, jpj 
     224               DO ji = 1, jpi 
     225                  zwz(ji,jj,mikt(ji,jj)) = 0.e0 
     226               END DO 
     227            END DO 
     228         ELSE 
     229            zwz(:,:,1) = 0.e0 
     230         END IF 
    217231         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions 
    218232         CALL lbc_lnk( zwz, 'W',  1. ) 
    … …  
    250264         END IF 
    251265         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    252          IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN   
    253            IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) ) + htr_adv(:) 
    254            IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) ) + str_adv(:) 
     266         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
     267           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) + htr_adv(:) 
     268           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) + str_adv(:) 
    255269         ENDIF 
    256270         ! 
    … …  
    358372 
    359373         ! upstream tracer flux in the k direction 
    360          ! Surface value 
    361          IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0._wp                        ! volume variable 
    362          ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface  
    363          ENDIF 
    364374         ! Interior value 
    365375         DO jk = 2, jpkm1 
    … …  
    372382            END DO 
    373383         END DO 
     384         ! Surface value 
     385         IF( lk_vvl ) THEN 
     386            IF ( ln_isfcav ) THEN 
     387               DO jj = 1, jpj 
     388                  DO ji = 1, jpi 
     389                     zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable +    isf 
     390                  END DO 
     391               END DO 
     392            ELSE 
     393               zwz(:,:,1) = 0.e0                              ! volume variable + no isf 
     394            END IF 
     395         ELSE 
     396            IF ( ln_isfcav ) THEN 
     397               DO jj = 1, jpj 
     398                  DO ji = 1, jpi 
     399                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface +    isf 
     400                  END DO 
     401               END DO 
     402            ELSE 
     403               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)                                               ! linear free surface + no isf 
     404            END IF 
     405         ENDIF 
    374406 
    375407         ! total advective trend 
    … …  
    398430         END IF 
    399431         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    400          IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN   
    401            IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) ) 
    402            IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) ) 
     432         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
     433           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) 
     434           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) 
    403435         ENDIF 
    404436 
    … …  
    524556         END IF 
    525557         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    526          IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN   
    527            IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) ) + htr_adv(:) 
    528            IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) ) + str_adv(:) 
     558         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
     559           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) + htr_adv(:) 
     560           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) + str_adv(:) 
    529561         ENDIF 
    530562         ! 
    … …  
    580612         &        paft * tmask + zbig * ( 1._wp - tmask )  ) 
    581613 
    582       DO jj = 2, jpjm1 
    583          DO ji = fs_2, fs_jpim1   ! vector opt. 
    584             DO jk = mikt(ji,jj), jpkm1 
    585                ikm1 = MAX(jk-1,mikt(ji,jj)) 
    586                z2dtt = p2dt(jk) 
    587                 
     614      DO jk = 1, jpkm1 
     615         ikm1 = MAX(jk-1,1) 
     616         z2dtt = p2dt(jk) 
     617         DO jj = 2, jpjm1 
     618            DO ji = fs_2, fs_jpim1   ! vector opt. 
     619 
    588620               ! search maximum in neighbourhood 
    589621               zup = MAX(  zbup(ji  ,jj  ,jk  ),   & 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90

    r4990 r5682  
    177177         END IF 
    178178         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes) 
    179          IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN   
    180             IF( jn == jp_tem )  htr_adv(:) = ptr_vj( ztv(:,:,:) ) 
    181             IF( jn == jp_sal )  str_adv(:) = ptr_vj( ztv(:,:,:) ) 
     179         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
     180            IF( jn == jp_tem )  htr_adv(:) = ptr_sj( ztv(:,:,:) ) 
     181            IF( jn == jp_sal )  str_adv(:) = ptr_sj( ztv(:,:,:) ) 
    182182         ENDIF 
    183183          

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/trabbc.F90

    r4990 r5682  
    2121   USE trdtra          ! trends manager: tracers  
    2222   USE in_out_manager  ! I/O manager 
     23   USE iom             ! I/O manager 
     24   USE fldread         ! read input fields 
     25   USE lbclnk            ! ocean lateral boundary conditions (or mpp link) 
     26   USE lib_mpp           ! distributed memory computing library 
    2327   USE prtctl          ! Print control 
    2428   USE wrk_nemo        ! Memory Allocation 
    … …  
    3741 
    3842   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   qgh_trd0   ! geothermal heating trend 
     43   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qgh              ! structure of input qgh (file informations, fields read) 
    3944  
    4045   !! * Substitutions 
    … …  
    4247   !!---------------------------------------------------------------------- 
    4348   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    44    !! $Id $  
     49   !! $Id$ 
    4550   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4651   !!---------------------------------------------------------------------- 
    … …  
    9297      END DO 
    9398      ! 
     99      CALL lbc_lnk( tsa(:,:,:,jp_tem) , 'T', 1. ) 
     100      ! 
    94101      IF( l_trdtra ) THEN        ! Save the geothermal heat flux trend for diagnostics 
    95102         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    … …  
    125132      INTEGER  ::   inum                ! temporary logical unit 
    126133      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    127       ! 
    128       NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst  
     134      INTEGER  ::   ierror              ! local integer 
     135      ! 
     136      TYPE(FLD_N)        ::   sn_qgh    ! informations about the geotherm. field to be read 
     137      CHARACTER(len=256) ::   cn_dir    ! Root directory for location of ssr files 
     138      ! 
     139      NAMELIST/nambbc/ln_trabbc, nn_geoflx, rn_geoflx_cst, sn_qgh, cn_dir  
    129140      !!---------------------------------------------------------------------- 
    130141 
    … …  
    161172         CASE ( 2 )                          !* variable geothermal heat flux : read the geothermal fluxes in mW/m2 
    162173            IF(lwp) WRITE(numout,*) '      *** variable geothermal heat flux' 
    163             CALL iom_open ( 'geothermal_heating.nc', inum ) 
    164             CALL iom_get  ( inum, jpdom_data, 'heatflow', qgh_trd0 ) 
    165             CALL iom_close( inum ) 
    166             qgh_trd0(:,:) = r1_rau0_rcp * qgh_trd0(:,:) * 1.e-3     ! conversion in W/m2 
     174            ! 
     175            ALLOCATE( sf_qgh(1), STAT=ierror ) 
     176            IF( ierror > 0 ) THEN 
     177               CALL ctl_stop( 'tra_bbc_init: unable to allocate sf_qgh structure' )   ; 
     178               RETURN 
     179            ENDIF 
     180            ALLOCATE( sf_qgh(1)%fnow(jpi,jpj,1)   ) 
     181            IF( sn_qgh%ln_tint )ALLOCATE( sf_qgh(1)%fdta(jpi,jpj,1,2) ) 
     182            ! fill sf_chl with sn_chl and control print 
     183            CALL fld_fill( sf_qgh, (/ sn_qgh /), cn_dir, 'tra_bbc_init',   & 
     184               &          'bottom temperature boundary condition', 'nambbc' ) 
     185 
     186            CALL fld_read( nit000, 1, sf_qgh )                         ! Read qgh data 
     187            qgh_trd0(:,:) = r1_rau0_rcp * sf_qgh(1)%fnow(:,:,1) * 1.e-3 ! conversion in W/m2 
    167188            ! 
    168189         CASE DEFAULT 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90

    r4990 r5682  
    2121   !!   tra_dmp       : update the tracer trend with the internal damping 
    2222   !!   tra_dmp_init  : initialization, namlist read, parameters control 
    23    !!   dtacof_zoom   : restoring coefficient for zoom domain 
    24    !!   dtacof        : restoring coefficient for global domain 
    25    !!   cofdis        : compute the distance to the coastline 
    2623   !!---------------------------------------------------------------------- 
    2724   USE oce            ! ocean: variables 
    … …  
    3936   USE wrk_nemo       ! Memory allocation 
    4037   USE timing         ! Timing 
     38   USE iom 
    4139 
    4240   IMPLICIT NONE 
    … …  
    4543   PUBLIC   tra_dmp      ! routine called by step.F90 
    4644   PUBLIC   tra_dmp_init ! routine called by opa.F90 
    47    PUBLIC   dtacof       ! routine called by tradmp.F90, trcdmp.F90 and dyndmp.F90 
    48    PUBLIC   dtacof_zoom  ! routine called by tradmp.F90, trcdmp.F90 and dyndmp.F90 
    49  
    50 !!gm  why all namelist variable public????   only ln_tradmp should be sufficient 
    5145 
    5246   !                               !!* Namelist namtra_dmp : T & S newtonian damping * 
     47   ! nn_zdmp and cn_resto are public as they are used by C1D/dyndmp.F90 
    5348   LOGICAL , PUBLIC ::   ln_tradmp   !: internal damping flag 
    54    INTEGER , PUBLIC ::   nn_hdmp     ! = 0/-1/'latitude' for damping over T and S 
    5549   INTEGER , PUBLIC ::   nn_zdmp     ! = 0/1/2 flag for damping in the mixed layer 
    56    REAL(wp), PUBLIC ::   rn_surf     ! surface time scale for internal damping        [days] 
    57    REAL(wp), PUBLIC ::   rn_bot      ! bottom time scale for internal damping         [days] 
    58    REAL(wp), PUBLIC ::   rn_dep      ! depth of transition between rn_surf and rn_bot [meters] 
    59    INTEGER , PUBLIC ::   nn_file     ! = 1 create a damping.coeff NetCDF file  
     50   CHARACTER(LEN=200) , PUBLIC :: cn_resto      ! name of netcdf file containing restoration coefficient field 
     51   ! 
     52 
    6053 
    6154   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   strdmp   !: damping salinity trend (psu/s) 
    … …  
    197190      !! ** Method  :   read the namtra_dmp namelist and check the parameters 
    198191      !!---------------------------------------------------------------------- 
    199       INTEGER  ::   ios   ! Local integer output status for namelist read 
    200       !! 
    201       NAMELIST/namtra_dmp/ ln_tradmp, nn_hdmp, nn_zdmp, rn_surf, rn_bot, rn_dep, nn_file 
    202       !!---------------------------------------------------------------------- 
    203       ! 
    204       REWIND( numnam_ref )              ! Namelist namtra_dmp in reference namelist : Temperature and salinity damping term 
     192      NAMELIST/namtra_dmp/ ln_tradmp, nn_zdmp, cn_resto 
     193      INTEGER ::  ios         ! Local integer for output status of namelist read 
     194      INTEGER :: imask        ! File handle  
     195      !! 
     196      !!---------------------------------------------------------------------- 
     197      ! 
     198      REWIND( numnam_ref )   ! Namelist namtra_dmp in reference namelist : T & S relaxation 
    205199      READ  ( numnam_ref, namtra_dmp, IOSTAT = ios, ERR = 901) 
    206200901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in reference namelist', lwp ) 
    207201      ! 
    208       REWIND( numnam_cfg )              ! Namelist namtra_dmp in configuration namelist : Temperature and salinity damping term 
     202      REWIND( numnam_cfg )   ! Namelist namtra_dmp in configuration namelist : T & S relaxation 
    209203      READ  ( numnam_cfg, namtra_dmp, IOSTAT = ios, ERR = 902 ) 
    210204902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist', lwp ) 
    211205      IF(lwm) WRITE ( numond, namtra_dmp ) 
    212        
    213       IF( lzoom .AND. .NOT. lk_c1d )   nn_zdmp = 0          ! restoring to climatology at closed north or south boundaries 
    214  
    215       IF(lwp) THEN                       ! Namelist print 
     206 
     207      IF(lwp) THEN                 !Namelist print 
    216208         WRITE(numout,*) 
    217          WRITE(numout,*) 'tra_dmp_init : T and S newtonian damping' 
     209         WRITE(numout,*) 'tra_dmp_init : T and S newtonian relaxation' 
    218210         WRITE(numout,*) '~~~~~~~' 
    219          WRITE(numout,*) '   Namelist namtra_dmp : set damping parameter' 
    220          WRITE(numout,*) '      add a damping term or not       ln_tradmp = ', ln_tradmp 
    221          WRITE(numout,*) '      T and S damping option          nn_hdmp   = ', nn_hdmp 
    222          WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp, '(non-C1D zoom: forced to 0)' 
    223          WRITE(numout,*) '      surface time scale (days)       rn_surf   = ', rn_surf 
    224          WRITE(numout,*) '      bottom time scale (days)        rn_bot    = ', rn_bot 
    225          WRITE(numout,*) '      depth of transition (meters)    rn_dep    = ', rn_dep 
    226          WRITE(numout,*) '      create a damping.coeff file     nn_file   = ', nn_file 
     211         WRITE(numout,*) '   Namelist namtra_dmp : set relaxation parameters' 
     212         WRITE(numout,*) '      Apply relaxation   or not       ln_tradmp = ', ln_tradmp 
     213         WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp 
     214         WRITE(numout,*) '      Damping file name               cn_resto  = ', cn_resto 
    227215         WRITE(numout,*) 
    228216      ENDIF 
    229217 
    230       IF( ln_tradmp ) THEN               ! initialization for T-S damping 
    231          ! 
     218      IF( ln_tradmp) THEN 
     219         ! 
     220         !Allocate arrays 
    232221         IF( tra_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_dmp_init: unable to allocate arrays' ) 
    233          ! 
    234 !!gm  I don't understand the specificities of c1d case......    
    235 !!gm  to be check with the autor of these lines 
    236           
    237 #if ! defined key_c1d 
    238          SELECT CASE ( nn_hdmp ) 
    239          CASE (  -1  )   ;   IF(lwp) WRITE(numout,*) '   tracer damping in the Med & Red seas only' 
    240          CASE ( 1:90 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping poleward of', nn_hdmp, ' degrees' 
    241          CASE DEFAULT 
    242             WRITE(ctmp1,*) '          bad flag value for nn_hdmp = ', nn_hdmp 
    243             CALL ctl_stop(ctmp1) 
     222 
     223         !Check values of nn_zdmp 
     224         SELECT CASE (nn_zdmp) 
     225         CASE ( 0 )  ; IF(lwp) WRITE(numout,*) '   tracer damping as specified by mask' 
     226         CASE ( 1 )  ; IF(lwp) WRITE(numout,*) '   no tracer damping in the turbocline' 
     227         CASE ( 2 )  ; IF(lwp) WRITE(numout,*) '   no tracer damping in the mixed layer' 
    244228         END SELECT 
    245          ! 
    246 #endif 
    247          SELECT CASE ( nn_zdmp ) 
    248          CASE ( 0 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping throughout the water column' 
    249          CASE ( 1 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the turbocline (avt > 5 cm2/s)' 
    250          CASE ( 2 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the mixed layer' 
    251          CASE DEFAULT 
    252             WRITE(ctmp1,*) 'bad flag value for nn_zdmp = ', nn_zdmp 
    253             CALL ctl_stop(ctmp1) 
    254          END SELECT 
    255          ! 
     229 
     230         !TG: Initialisation of dtatsd - Would it be better to have dmpdta routine 
     231         !so can damp to something other than intitial conditions files? 
    256232         IF( .NOT.ln_tsd_tradmp ) THEN 
    257233            CALL ctl_warn( 'tra_dmp_init: read T-S data not initialized, we force ln_tsd_tradmp=T' ) 
    258234            CALL dta_tsd_init( ld_tradmp=ln_tradmp )        ! forces the initialisation of T-S data 
    259235         ENDIF 
    260          ! 
    261          strdmp(:,:,:) = 0._wp       ! internal damping salinity trend (used in asmtrj) 
     236 
     237         !initialise arrays - Are these actually used anywhere else? 
     238         strdmp(:,:,:) = 0._wp 
    262239         ttrdmp(:,:,:) = 0._wp 
    263          !                          ! Damping coefficients initialization 
    264          IF( lzoom .AND. .NOT. lk_c1d ) THEN   ;   CALL dtacof_zoom( resto ) 
    265          ELSE               ;   CALL dtacof( nn_hdmp, rn_surf, rn_bot, rn_dep, nn_file, 'TRA', resto ) 
    266          ENDIF 
    267          ! 
    268       ENDIF 
    269       ! 
     240 
     241         !Read in mask from file 
     242         CALL iom_open ( cn_resto, imask) 
     243         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto) 
     244         CALL iom_close( imask ) 
     245       ENDIF 
     246 
    270247   END SUBROUTINE tra_dmp_init 
    271248 
    272  
    273    SUBROUTINE dtacof_zoom( presto ) 
    274       !!---------------------------------------------------------------------- 
    275       !!                  ***  ROUTINE dtacof_zoom  *** 
    276       !! 
    277       !! ** Purpose :   Compute the damping coefficient for zoom domain 
    278       !! 
    279       !! ** Method  : - set along closed boundary due to zoom a damping over 
    280       !!                6 points with a max time scale of 5 days. 
    281       !!              - ORCA arctic/antarctic zoom: set the damping along 
    282       !!                south/north boundary over a latitude strip. 
    283       !! 
    284       !! ** Action  : - resto, the damping coeff. for T and S 
    285       !!---------------------------------------------------------------------- 
    286       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)  ::   presto   ! restoring coeff. (s-1) 
    287       ! 
    288       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    289       REAL(wp) ::   zlat, zlat0, zlat1, zlat2, z1_5d   ! local scalar 
    290       REAL(wp), DIMENSION(6)  ::   zfact               ! 1Dworkspace 
    291       !!---------------------------------------------------------------------- 
    292       ! 
    293       IF( nn_timing == 1 )  CALL timing_start( 'dtacof_zoom') 
    294       ! 
    295  
    296       zfact(1) =  1._wp 
    297       zfact(2) =  1._wp 
    298       zfact(3) = 11._wp / 12._wp 
    299       zfact(4) =  8._wp / 12._wp 
    300       zfact(5) =  4._wp / 12._wp 
    301       zfact(6) =  1._wp / 12._wp 
    302       zfact(:) = zfact(:) / ( 5._wp * rday )    ! 5 days max restoring time scale 
    303  
    304       presto(:,:,:) = 0._wp 
    305  
    306       ! damping along the forced closed boundary over 6 grid-points 
    307       DO jn = 1, 6 
    308          IF( lzoom_w )   presto( mi0(jn+jpizoom):mi1(jn+jpizoom), : , : )                    = zfact(jn)   ! west  closed 
    309          IF( lzoom_s )   presto( : , mj0(jn+jpjzoom):mj1(jn+jpjzoom), : )                    = zfact(jn)   ! south closed  
    310          IF( lzoom_e )   presto( mi0(jpiglo+jpizoom-1-jn):mi1(jpiglo+jpizoom-1-jn) , : , : ) = zfact(jn)   ! east  closed  
    311          IF( lzoom_n )   presto( : , mj0(jpjglo+jpjzoom-1-jn):mj1(jpjglo+jpjzoom-1-jn) , : ) = zfact(jn)   ! north closed 
    312       END DO 
    313  
    314       !                                           ! ==================================================== 
    315       IF( cp_cfz == "arctic" .OR. cp_cfz == "antarctic" ) THEN   !  ORCA configuration : arctic or antarctic zoom 
    316          !                                        ! ==================================================== 
    317          IF(lwp) WRITE(numout,*) 
    318          IF(lwp .AND. cp_cfz == "arctic" ) WRITE(numout,*) '              dtacof_zoom : ORCA    Arctic zoom' 
    319          IF(lwp .AND. cp_cfz == "antarctic" ) WRITE(numout,*) '           dtacof_zoom : ORCA Antarctic zoom' 
    320          IF(lwp) WRITE(numout,*) 
    321          ! 
    322          !                          ! Initialization :  
    323          presto(:,:,:) = 0._wp 
    324          zlat0 = 10._wp                     ! zlat0 : latitude strip where resto decreases 
    325          zlat1 = 30._wp                     ! zlat1 : resto = 1 before zlat1 
    326          zlat2 = zlat1 + zlat0              ! zlat2 : resto decreases from 1 to 0 between zlat1 and zlat2 
    327          z1_5d = 1._wp / ( 5._wp * rday )   ! z1_5d : 1 / 5days 
    328  
    329          DO jk = 2, jpkm1           ! Compute arrays resto ; value for internal damping : 5 days 
    330             DO jj = 1, jpj 
    331                DO ji = 1, jpi 
    332                   zlat = ABS( gphit(ji,jj) ) 
    333                   IF( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN 
    334                      presto(ji,jj,jk) = 0.5_wp * z1_5d * (  1._wp - COS( rpi*(zlat2-zlat)/zlat0 )  )  
    335                   ELSEIF( zlat < zlat1 ) THEN 
    336                      presto(ji,jj,jk) = z1_5d 
    337                   ENDIF 
    338                END DO 
    339             END DO 
    340          END DO 
    341          ! 
    342       ENDIF 
    343       !                             ! Mask resto array 
    344       presto(:,:,:) = presto(:,:,:) * tmask(:,:,:) 
    345       ! 
    346       IF( nn_timing == 1 )  CALL timing_stop( 'dtacof_zoom') 
    347       ! 
    348    END SUBROUTINE dtacof_zoom 
    349  
    350  
    351    SUBROUTINE dtacof( kn_hdmp, pn_surf, pn_bot, pn_dep,  & 
    352       &               kn_file, cdtype , presto           ) 
    353       !!---------------------------------------------------------------------- 
    354       !!                  ***  ROUTINE dtacof  *** 
    355       !! 
    356       !! ** Purpose :   Compute the damping coefficient 
    357       !! 
    358       !! ** Method  :   Arrays defining the damping are computed for each grid 
    359       !!                point for temperature and salinity (resto) 
    360       !!                Damping depends on distance to coast, depth and latitude 
    361       !! 
    362       !! ** Action  : - resto, the damping coeff. for T and S 
    363       !!---------------------------------------------------------------------- 
    364       USE iom 
    365       USE ioipsl 
    366       !! 
    367       INTEGER                         , INTENT(in   )  ::  kn_hdmp    ! damping option 
    368       REAL(wp)                        , INTENT(in   )  ::  pn_surf    ! surface time scale (days) 
    369       REAL(wp)                        , INTENT(in   )  ::  pn_bot     ! bottom time scale (days) 
    370       REAL(wp)                        , INTENT(in   )  ::  pn_dep     ! depth of transition (meters) 
    371       INTEGER                         , INTENT(in   )  ::  kn_file    ! save the damping coef on a file or not 
    372       CHARACTER(len=3)                , INTENT(in   )  ::  cdtype     ! =TRA, TRC or DYN (tracer/dynamics indicator) 
    373       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)  ::  presto     ! restoring coeff. (s-1) 
    374       ! 
    375       INTEGER  ::   ji, jj, jk                  ! dummy loop indices 
    376       INTEGER  ::   ii0, ii1, ij0, ij1          ! local integers 
    377       INTEGER  ::   inum0, icot                 !   -       - 
    378       REAL(wp) ::   zinfl, zlon                 ! local scalars 
    379       REAL(wp) ::   zlat, zlat0, zlat1, zlat2   !   -      - 
    380       REAL(wp) ::   zsdmp, zbdmp                !   -      - 
    381       CHARACTER(len=20)                   :: cfile 
    382       REAL(wp), POINTER, DIMENSION(:    ) :: zhfac  
    383       REAL(wp), POINTER, DIMENSION(:,:  ) :: zmrs  
    384       REAL(wp), POINTER, DIMENSION(:,:,:) :: zdct  
    385       !!---------------------------------------------------------------------- 
    386       ! 
    387       IF( nn_timing == 1 )  CALL timing_start('dtacof') 
    388       ! 
    389       CALL wrk_alloc( jpk, zhfac          ) 
    390       CALL wrk_alloc( jpi, jpj, zmrs      ) 
    391       CALL wrk_alloc( jpi, jpj, jpk, zdct ) 
    392 #if defined key_c1d 
    393       !                                   ! ==================== 
    394       !                                   !  C1D configuration : local domain 
    395       !                                   ! ==================== 
    396       ! 
    397       IF(lwp) WRITE(numout,*) 
    398       IF(lwp) WRITE(numout,*) '              dtacof : C1D 3x3 local domain' 
    399       IF(lwp) WRITE(numout,*) '              -----------------------------' 
    400       ! 
    401       presto(:,:,:) = 0._wp 
    402       ! 
    403       zsdmp = 1._wp / ( pn_surf * rday ) 
    404       zbdmp = 1._wp / ( pn_bot  * rday ) 
    405       DO jk = 2, jpkm1 
    406          DO jj = 1, jpj 
    407             DO ji = 1, jpi 
    408                !   ONLY vertical variation from zsdmp (sea surface) to zbdmp (bottom) 
    409                presto(ji,jj,jk) = zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(ji,jj,jk)/pn_dep) 
    410             END DO 
    411          END DO 
    412       END DO 
    413       ! 
    414       presto(:,:, : ) = presto(:,:,:) * tmask(:,:,:) 
    415 #else 
    416       !                                   ! ==================== 
    417       !                                   !  ORCA configuration : global domain 
    418       !                                   ! ==================== 
    419       ! 
    420       IF(lwp) WRITE(numout,*) 
    421       IF(lwp) WRITE(numout,*) '              dtacof : Global domain of ORCA' 
    422       IF(lwp) WRITE(numout,*) '              ------------------------------' 
    423       ! 
    424       presto(:,:,:) = 0._wp 
    425       ! 
    426       IF( kn_hdmp > 0 ) THEN      !  Damping poleward of 'nn_hdmp' degrees  ! 
    427          !                        !-----------------------------------------! 
    428          IF(lwp) WRITE(numout,*) 
    429          IF(lwp) WRITE(numout,*) '              Damping poleward of ', kn_hdmp, ' deg.' 
    430          ! 
    431          CALL iom_open ( 'dist.coast.nc', icot, ldstop = .FALSE. ) 
    432          ! 
    433          IF( icot > 0 ) THEN          ! distance-to-coast read in file 
    434             CALL iom_get  ( icot, jpdom_data, 'Tcoast', zdct ) 
    435             CALL iom_close( icot ) 
    436          ELSE                         ! distance-to-coast computed and saved in file (output in zdct) 
    437             CALL cofdis( zdct ) 
    438          ENDIF 
    439  
    440          !                            ! Compute arrays resto  
    441          zinfl = 1000.e3_wp                ! distance of influence for damping term 
    442          zlat0 = 10._wp                    ! latitude strip where resto decreases 
    443          zlat1 = REAL( kn_hdmp )           ! resto = 0 between -zlat1 and zlat1 
    444          zlat2 = zlat1 + zlat0             ! resto increases from 0 to 1 between |zlat1| and |zlat2| 
    445  
    446          DO jj = 1, jpj 
    447             DO ji = 1, jpi 
    448                zlat = ABS( gphit(ji,jj) ) 
    449                IF ( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN 
    450                   presto(ji,jj,1) = 0.5_wp * (  1._wp - COS( rpi*(zlat-zlat1)/zlat0 )  ) 
    451                ELSEIF ( zlat > zlat2 ) THEN 
    452                   presto(ji,jj,1) = 1._wp 
    453                ENDIF 
    454             END DO 
    455          END DO 
    456  
    457          IF ( kn_hdmp == 20 ) THEN       ! North Indian ocean (20N/30N x 45E/100E) : resto=0 
    458             DO jj = 1, jpj 
    459                DO ji = 1, jpi 
    460                   zlat = gphit(ji,jj) 
    461                   zlon = MOD( glamt(ji,jj), 360._wp ) 
    462                   IF ( zlat1 < zlat .AND. zlat < zlat2 .AND. 45._wp < zlon .AND. zlon < 100._wp ) THEN 
    463                      presto(ji,jj,1) = 0._wp 
    464                   ENDIF 
    465                END DO 
    466             END DO 
    467          ENDIF 
    468  
    469          zsdmp = 1._wp / ( pn_surf * rday ) 
    470          zbdmp = 1._wp / ( pn_bot  * rday ) 
    471          DO jk = 2, jpkm1 
    472             DO jj = 1, jpj 
    473                DO ji = 1, jpi 
    474                   zdct(ji,jj,jk) = MIN( zinfl, zdct(ji,jj,jk) ) 
    475                   !   ... Decrease the value in the vicinity of the coast 
    476                   presto(ji,jj,jk) = presto(ji,jj,1 ) * 0.5_wp * (  1._wp - COS( rpi*zdct(ji,jj,jk)/zinfl)  ) 
    477                   !   ... Vertical variation from zsdmp (sea surface) to zbdmp (bottom) 
    478                   presto(ji,jj,jk) = presto(ji,jj,jk) * (  zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(ji,jj,jk)/pn_dep)  ) 
    479                END DO 
    480             END DO 
    481          END DO 
    482          ! 
    483       ENDIF 
    484  
    485       !                                  ! ========================= 
    486       !                                  !  Med and Red Sea damping    (ORCA configuration only) 
    487       !                                  ! ========================= 
    488       IF( cp_cfg == "orca" .AND. ( kn_hdmp > 0 .OR. kn_hdmp == -1 ) ) THEN 
    489          IF(lwp)WRITE(numout,*) 
    490          IF(lwp)WRITE(numout,*) '              ORCA configuration: Damping in Med and Red Seas' 
    491          ! 
    492          zmrs(:,:) = 0._wp 
    493          ! 
    494          SELECT CASE ( jp_cfg ) 
    495          !                                           ! ======================= 
    496          CASE ( 4 )                                  !  ORCA_R4 configuration  
    497             !                                        ! ======================= 
    498             ij0 =  50   ;   ij1 =  56                    ! Mediterranean Sea 
    499  
    500             ii0 =  81   ;   ii1 =  91   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    501             ij0 =  50   ;   ij1 =  55 
    502             ii0 =  75   ;   ii1 =  80   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    503             ij0 =  52   ;   ij1 =  53 
    504             ii0 =  70   ;   ii1 =  74   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    505             ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea 
    506             DO jk = 1, 17 
    507                zhfac (jk) = 0.5_wp * (  1._wp - COS( rpi * REAL(jk-1,wp) / 16._wp )  ) / rday 
    508             END DO 
    509             DO jk = 18, jpkm1 
    510                zhfac (jk) = 1._wp / rday 
    511             END DO 
    512             !                                        ! ======================= 
    513          CASE ( 2 )                                  !  ORCA_R2 configuration  
    514             !                                        ! ======================= 
    515             ij0 =  96   ;   ij1 = 110                    ! Mediterranean Sea 
    516             ii0 = 157   ;   ii1 = 181   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    517             ij0 = 100   ;   ij1 = 110 
    518             ii0 = 144   ;   ii1 = 156   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    519             ij0 = 100   ;   ij1 = 103 
    520             ii0 = 139   ;   ii1 = 143   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    521             ! 
    522             ij0 = 101   ;   ij1 = 102                    ! Decrease before Gibraltar Strait 
    523             ii0 = 139   ;   ii1 = 141   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0._wp 
    524             ii0 = 142   ;   ii1 = 142   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp 
    525             ii0 = 143   ;   ii1 = 143   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp 
    526             ii0 = 144   ;   ii1 = 144   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75_wp 
    527             ! 
    528             ij0 =  87   ;   ij1 =  96                    ! Red Sea 
    529             ii0 = 147   ;   ii1 = 163   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    530             ! 
    531             ij0 =  91   ;   ij1 =  91                    ! Decrease before Bab el Mandeb Strait 
    532             ii0 = 153   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.80_wp 
    533             ij0 =  90   ;   ij1 =  90 
    534             ii0 = 153   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp 
    535             ij0 =  89   ;   ij1 =  89 
    536             ii0 = 158   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp 
    537             ij0 =  88   ;   ij1 =  88 
    538             ii0 = 160   ;   ii1 = 163   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0._wp 
    539             ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea 
    540             DO jk = 1, 17 
    541                zhfac (jk) = 0.5_wp * (  1._wp - COS( rpi * REAL(jk-1,wp) / 16._wp )  ) / rday 
    542             END DO 
    543             DO jk = 18, jpkm1 
    544                zhfac (jk) = 1._wp / rday 
    545             END DO 
    546             !                                        ! ======================= 
    547          CASE ( 05 )                                 !  ORCA_R05 configuration 
    548             !                                        ! ======================= 
    549             ii0 = 568   ;   ii1 = 574                    ! Mediterranean Sea 
    550             ij0 = 324   ;   ij1 = 333   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    551             ii0 = 575   ;   ii1 = 658 
    552             ij0 = 314   ;   ij1 = 366   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    553             ! 
    554             ii0 = 641   ;   ii1 = 651                    ! Black Sea (remaining part 
    555             ij0 = 367   ;   ij1 = 372   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    556             ! 
    557             ij0 = 324   ;   ij1 = 333                    ! Decrease before Gibraltar Strait 
    558             ii0 = 565   ;   ii1 = 565   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp 
    559             ii0 = 566   ;   ii1 = 566   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp 
    560             ii0 = 567   ;   ii1 = 567   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75_wp 
    561             ! 
    562             ii0 = 641   ;   ii1 = 665                    ! Red Sea 
    563             ij0 = 270   ;   ij1 = 310   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    564             ! 
    565             ii0 = 666   ;   ii1 = 675                    ! Decrease before Bab el Mandeb Strait 
    566             ij0 = 270   ;   ij1 = 290    
    567             DO ji = mi0(ii0), mi1(ii1) 
    568                zmrs( ji , mj0(ij0):mj1(ij1) ) = 0.1_wp * ABS( FLOAT(ji - mi1(ii1)) ) 
    569             END DO  
    570             zsdmp = 1._wp / ( pn_surf * rday ) 
    571             zbdmp = 1._wp / ( pn_bot  * rday ) 
    572             DO jk = 1, jpk 
    573                zhfac(jk) = (  zbdmp + (zsdmp-zbdmp) * EXP( -fsdept(1,1,jk)/pn_dep )  ) 
    574             END DO 
    575             !                                       ! ======================== 
    576          CASE ( 025 )                               !  ORCA_R025 configuration  
    577             !                                       ! ======================== 
    578             CALL ctl_stop( ' Not yet implemented in ORCA_R025' ) 
    579             ! 
    580          END SELECT 
    581  
    582          DO jk = 1, jpkm1 
    583             presto(:,:,jk) = zmrs(:,:) * zhfac(jk) + ( 1._wp - zmrs(:,:) ) * presto(:,:,jk) 
    584          END DO 
    585  
    586          ! Mask resto array and set to 0 first and last levels 
    587          presto(:,:, : ) = presto(:,:,:) * tmask(:,:,:) 
    588          presto(:,:, 1 ) = 0._wp 
    589          presto(:,:,jpk) = 0._wp 
    590          !                         !--------------------! 
    591       ELSE                         !     No damping     ! 
    592          !                         !--------------------! 
    593          CALL ctl_stop( 'Choose a correct value of nn_hdmp or put ln_tradmp to FALSE' ) 
    594       ENDIF 
    595 #endif 
    596  
    597       !                            !--------------------------------! 
    598       IF( kn_file == 1 ) THEN      !  save damping coef. in a file  ! 
    599          !                         !--------------------------------! 
    600          IF(lwp) WRITE(numout,*) '              create damping.coeff.nc file' 
    601          IF( cdtype == 'TRA' ) cfile = 'damping.coeff' 
    602          IF( cdtype == 'TRC' ) cfile = 'damping.coeff.trc' 
    603          IF( cdtype == 'DYN' ) cfile = 'damping.coeff.dyn' 
    604          cfile = TRIM( cfile ) 
    605          CALL iom_open  ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib ) 
    606          CALL iom_rstput( 0, 0, inum0, 'Resto', presto ) 
    607          CALL iom_close ( inum0 ) 
    608       ENDIF 
    609       ! 
    610       CALL wrk_dealloc( jpk, zhfac) 
    611       CALL wrk_dealloc( jpi, jpj, zmrs ) 
    612       CALL wrk_dealloc( jpi, jpj, jpk, zdct ) 
    613       ! 
    614       IF( nn_timing == 1 )  CALL timing_stop('dtacof') 
    615       ! 
    616    END SUBROUTINE dtacof 
    617  
    618  
    619    SUBROUTINE cofdis( pdct ) 
    620       !!---------------------------------------------------------------------- 
    621       !!                 ***  ROUTINE cofdis  *** 
    622       !! 
    623       !! ** Purpose :   Compute the distance between ocean T-points and the 
    624       !!      ocean model coastlines. Save the distance in a NetCDF file. 
    625       !! 
    626       !! ** Method  :   For each model level, the distance-to-coast is  
    627       !!      computed as follows :  
    628       !!       - The coastline is defined as the serie of U-,V-,F-points 
    629       !!      that are at the ocean-land bound. 
    630       !!       - For each ocean T-point, the distance-to-coast is then  
    631       !!      computed as the smallest distance (on the sphere) between the  
    632       !!      T-point and all the coastline points. 
    633       !!       - For land T-points, the distance-to-coast is set to zero. 
    634       !!      C A U T I O N : Computation not yet implemented in mpp case. 
    635       !! 
    636       !! ** Action  : - pdct, distance to the coastline (argument) 
    637       !!              - NetCDF file 'dist.coast.nc'  
    638       !!---------------------------------------------------------------------- 
    639       USE ioipsl      ! IOipsl librairy 
    640       !! 
    641       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) ::   pdct   ! distance to the coastline 
    642       !! 
    643       INTEGER ::   ji, jj, jk, jl   ! dummy loop indices 
    644       INTEGER ::   iju, ijt, icoast, itime, ierr, icot   ! local integers 
    645       CHARACTER (len=32) ::   clname                     ! local name 
    646       REAL(wp) ::   zdate0                               ! local scalar 
    647       REAL(wp), POINTER, DIMENSION(:,:) ::  zxt, zyt, zzt, zmask 
    648       REAL(wp), POINTER, DIMENSION(:  ) ::  zxc, zyc, zzc, zdis    ! temporary workspace 
    649       LOGICAL , ALLOCATABLE, DIMENSION(:,:) ::  llcotu, llcotv, llcotf   ! 2D logical workspace 
    650       !!---------------------------------------------------------------------- 
    651       ! 
    652       IF( nn_timing == 1 )  CALL timing_start('cofdis') 
    653       ! 
    654       CALL wrk_alloc( jpi, jpj , zxt, zyt, zzt, zmask    ) 
    655       CALL wrk_alloc( 3*jpi*jpj, zxc, zyc, zzc, zdis     ) 
    656       ALLOCATE( llcotu(jpi,jpj), llcotv(jpi,jpj), llcotf(jpi,jpj)  ) 
    657       ! 
    658       IF( lk_mpp    )   CALL mpp_sum( ierr ) 
    659       IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'cofdis: requested local arrays unavailable') 
    660  
    661       ! 0. Initialization 
    662       ! ----------------- 
    663       IF(lwp) WRITE(numout,*) 
    664       IF(lwp) WRITE(numout,*) 'cofdis : compute the distance to coastline' 
    665       IF(lwp) WRITE(numout,*) '~~~~~~' 
    666       IF(lwp) WRITE(numout,*) 
    667       IF( lk_mpp ) & 
    668            & CALL ctl_stop('         Computation not yet implemented with key_mpp_...', & 
    669            &               '         Rerun the code on another computer or ', & 
    670            &               '         create the "dist.coast.nc" file using IDL' ) 
    671  
    672       pdct(:,:,:) = 0._wp 
    673       zxt(:,:) = COS( rad * gphit(:,:) ) * COS( rad * glamt(:,:) ) 
    674       zyt(:,:) = COS( rad * gphit(:,:) ) * SIN( rad * glamt(:,:) ) 
    675       zzt(:,:) = SIN( rad * gphit(:,:) ) 
    676  
    677  
    678       ! 1. Loop on vertical levels 
    679       ! -------------------------- 
    680       !                                                ! =============== 
    681       DO jk = 1, jpkm1                                 ! Horizontal slab 
    682          !                                             ! =============== 
    683          ! Define the coastline points (U, V and F) 
    684          DO jj = 2, jpjm1 
    685             DO ji = 2, jpim1 
    686                zmask(ji,jj) =  ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 
    687                    &           + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) ) 
    688                llcotu(ji,jj) = ( tmask(ji,jj,  jk) + tmask(ji+1,jj  ,jk) == 1._wp )  
    689                llcotv(ji,jj) = ( tmask(ji,jj  ,jk) + tmask(ji  ,jj+1,jk) == 1._wp )  
    690                llcotf(ji,jj) = ( zmask(ji,jj) > 0._wp ) .AND. ( zmask(ji,jj) < 4._wp ) 
    691             END DO 
    692          END DO 
    693  
    694          ! Lateral boundaries conditions 
    695          llcotu(:, 1 ) = umask(:,  2  ,jk) == 1 
    696          llcotu(:,jpj) = umask(:,jpjm1,jk) == 1 
    697          llcotv(:, 1 ) = vmask(:,  2  ,jk) == 1 
    698          llcotv(:,jpj) = vmask(:,jpjm1,jk) == 1 
    699          llcotf(:, 1 ) = fmask(:,  2  ,jk) == 1 
    700          llcotf(:,jpj) = fmask(:,jpjm1,jk) == 1 
    701  
    702          IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 
    703             llcotu( 1 ,:) = llcotu(jpim1,:) 
    704             llcotu(jpi,:) = llcotu(  2  ,:) 
    705             llcotv( 1 ,:) = llcotv(jpim1,:) 
    706             llcotv(jpi,:) = llcotv(  2  ,:) 
    707             llcotf( 1 ,:) = llcotf(jpim1,:) 
    708             llcotf(jpi,:) = llcotf(  2  ,:) 
    709          ELSE 
    710             llcotu( 1 ,:) = umask(  2  ,:,jk) == 1 
    711             llcotu(jpi,:) = umask(jpim1,:,jk) == 1 
    712             llcotv( 1 ,:) = vmask(  2  ,:,jk) == 1 
    713             llcotv(jpi,:) = vmask(jpim1,:,jk) == 1 
    714             llcotf( 1 ,:) = fmask(  2  ,:,jk) == 1 
    715             llcotf(jpi,:) = fmask(jpim1,:,jk) == 1 
    716          ENDIF 
    717          IF( nperio == 3 .OR. nperio == 4 ) THEN 
    718             DO ji = 1, jpim1 
    719                iju = jpi - ji + 1 
    720                llcotu(ji,jpj  ) = llcotu(iju,jpj-2) 
    721                llcotf(ji,jpjm1) = llcotf(iju,jpj-2) 
    722                llcotf(ji,jpj  ) = llcotf(iju,jpj-3) 
    723             END DO 
    724             DO ji = jpi/2, jpim1 
    725                iju = jpi - ji + 1 
    726                llcotu(ji,jpjm1) = llcotu(iju,jpjm1) 
    727             END DO 
    728             DO ji = 2, jpi 
    729                ijt = jpi - ji + 2 
    730                llcotv(ji,jpjm1) = llcotv(ijt,jpj-2) 
    731                llcotv(ji,jpj  ) = llcotv(ijt,jpj-3) 
    732             END DO 
    733          ENDIF 
    734          IF( nperio == 5 .OR. nperio == 6 ) THEN 
    735             DO ji = 1, jpim1 
    736                iju = jpi - ji 
    737                llcotu(ji,jpj  ) = llcotu(iju,jpjm1) 
    738                llcotf(ji,jpj  ) = llcotf(iju,jpj-2) 
    739             END DO 
    740             DO ji = jpi/2, jpim1 
    741                iju = jpi - ji 
    742                llcotf(ji,jpjm1) = llcotf(iju,jpjm1) 
    743             END DO 
    744             DO ji = 1, jpi 
    745                ijt = jpi - ji + 1 
    746                llcotv(ji,jpj  ) = llcotv(ijt,jpjm1) 
    747             END DO 
    748             DO ji = jpi/2+1, jpi 
    749                ijt = jpi - ji + 1 
    750                llcotv(ji,jpjm1) = llcotv(ijt,jpjm1) 
    751             END DO 
    752          ENDIF 
    753  
    754          ! Compute cartesian coordinates of coastline points 
    755          ! and the number of coastline points 
    756          icoast = 0 
    757          DO jj = 1, jpj 
    758             DO ji = 1, jpi 
    759                IF( llcotf(ji,jj) ) THEN 
    760                   icoast = icoast + 1 
    761                   zxc(icoast) = COS( rad*gphif(ji,jj) ) * COS( rad*glamf(ji,jj) ) 
    762                   zyc(icoast) = COS( rad*gphif(ji,jj) ) * SIN( rad*glamf(ji,jj) ) 
    763                   zzc(icoast) = SIN( rad*gphif(ji,jj) ) 
    764                ENDIF 
    765                IF( llcotu(ji,jj) ) THEN 
    766                   icoast = icoast+1 
    767                   zxc(icoast) = COS( rad*gphiu(ji,jj) ) * COS( rad*glamu(ji,jj) ) 
    768                   zyc(icoast) = COS( rad*gphiu(ji,jj) ) * SIN( rad*glamu(ji,jj) ) 
    769                   zzc(icoast) = SIN( rad*gphiu(ji,jj) ) 
    770                ENDIF 
    771                IF( llcotv(ji,jj) ) THEN 
    772                   icoast = icoast+1 
    773                   zxc(icoast) = COS( rad*gphiv(ji,jj) ) * COS( rad*glamv(ji,jj) ) 
    774                   zyc(icoast) = COS( rad*gphiv(ji,jj) ) * SIN( rad*glamv(ji,jj) ) 
    775                   zzc(icoast) = SIN( rad*gphiv(ji,jj) ) 
    776                ENDIF 
    777             END DO 
    778          END DO 
    779  
    780          ! Distance for the T-points 
    781          DO jj = 1, jpj 
    782             DO ji = 1, jpi 
    783                IF( tmask(ji,jj,jk) == 0._wp ) THEN 
    784                   pdct(ji,jj,jk) = 0._wp 
    785                ELSE 
    786                   DO jl = 1, icoast 
    787                      zdis(jl) = ( zxt(ji,jj) - zxc(jl) )**2   & 
    788                         &     + ( zyt(ji,jj) - zyc(jl) )**2   & 
    789                         &     + ( zzt(ji,jj) - zzc(jl) )**2 
    790                   END DO 
    791                   pdct(ji,jj,jk) = ra * SQRT( MINVAL( zdis(1:icoast) ) ) 
    792                ENDIF 
    793             END DO 
    794          END DO 
    795          !                                                ! =============== 
    796       END DO                                              !   End of slab 
    797       !                                                   ! =============== 
    798  
    799  
    800       ! 2. Create the  distance to the coast file in NetCDF format 
    801       ! ----------------------------------------------------------     
    802       clname = 'dist.coast' 
    803       itime  = 0 
    804       CALL ymds2ju( 0     , 1       , 1     , 0._wp , zdate0 ) 
    805       CALL restini( 'NONE', jpi     , jpj   , glamt, gphit ,   & 
    806          &          jpk   , gdept_1d, clname, itime, zdate0,   & 
    807          &          rdt   , icot                         ) 
    808       CALL restput( icot, 'Tcoast', jpi, jpj, jpk, 0, pdct ) 
    809       CALL restclo( icot ) 
    810       ! 
    811       CALL wrk_dealloc( jpi, jpj , zxt, zyt, zzt, zmask    ) 
    812       CALL wrk_dealloc( 3*jpi*jpj, zxc, zyc, zzc, zdis     ) 
    813       DEALLOCATE( llcotu, llcotv, llcotf  ) 
    814       ! 
    815       IF( nn_timing == 1 )  CALL timing_stop('cofdis') 
    816       ! 
    817    END SUBROUTINE cofdis 
    818    !!====================================================================== 
    819249END MODULE tradmp 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

    r4990 r5682  
    290290      IF(lwp) WRITE(numout,*) '              homogeneous ocean T = ', zt0, ' S = ',zs0 
    291291 
     292      ! Initialisation of gtui/gtvi in case of no cavity 
     293      IF ( .NOT. ln_isfcav ) THEN 
     294         gtui(:,:,:) = 0.0_wp 
     295         gtvi(:,:,:) = 0.0_wp 
     296      END IF 
    292297      !                                        ! T & S profile (to be coded +namelist parameter 
    293298 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90

    r4990 r5682  
    116116               END DO 
    117117            END DO 
    118  
    119118            !                          !==  Laplacian  ==! 
    120119            ! 
    … …  
    125124               END DO 
    126125            END DO 
     126            ! 
    127127            IF( ln_zps ) THEN                ! set gradient at partial step level (last ocean level) 
    128128               DO jj = 1, jpjm1 
    … …  
    130130                     IF( mbku(ji,jj) == jk )  ztu(ji,jj,jk) = zeeu(ji,jj) * pgu(ji,jj,jn) 
    131131                     IF( mbkv(ji,jj) == jk )  ztv(ji,jj,jk) = zeev(ji,jj) * pgv(ji,jj,jn) 
    132                      ! (ISH) 
    133                      IF( miku(ji,jj) == jk )  ztu(ji,jj,jk) = zeeu(ji,jj) * pgui(ji,jj,jn) 
    134                      IF( mikv(ji,jj) == jk )  ztu(ji,jj,jk) = zeev(ji,jj) * pgvi(ji,jj,jn) 
    135132                  END DO 
    136133               END DO 
    137134            ENDIF 
     135            ! (ISH) 
     136            IF( ln_zps .AND. ln_isfcav ) THEN ! set gradient at partial step level (first ocean level in a cavity) 
     137               DO jj = 1, jpjm1 
     138                  DO ji = 1, jpim1 
     139                     IF( miku(ji,jj) == MAX(jk,2) )  ztu(ji,jj,jk) = zeeu(ji,jj) * pgui(ji,jj,jn) 
     140                     IF( mikv(ji,jj) == MAX(jk,2) )  ztu(ji,jj,jk) = zeev(ji,jj) * pgvi(ji,jj,jn) 
     141                  END DO 
     142               END DO 
     143            ENDIF 
     144            ! 
    138145            DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient 
    139146               DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    166173         !                                                 
    167174         ! "zonal" mean lateral diffusive heat and salt transport 
    168          IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN   
    169            IF( jn == jp_tem )  htr_ldf(:) = ptr_vj( ztv(:,:,:) ) 
    170            IF( jn == jp_sal )  str_ldf(:) = ptr_vj( ztv(:,:,:) ) 
     175         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN   
     176           IF( jn == jp_tem )  htr_ldf(:) = ptr_sj( ztv(:,:,:) ) 
     177           IF( jn == jp_sal )  str_ldf(:) = ptr_sj( ztv(:,:,:) ) 
    171178         ENDIF 
    172179         !                                                ! =========== 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90

    r4292 r5682  
    247247         !                                                ! =============== 
    248248         ! "Poleward" diffusive heat or salt transport 
    249          IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( kaht == 2 ) .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
     249         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( kaht == 2 ) ) THEN 
    250250            ! note sign is reversed to give down-gradient diffusive transports (#1043) 
    251             IF( jn == jp_tem)   htr_ldf(:) = ptr_vj( -zftv(:,:,:) ) 
    252             IF( jn == jp_sal)   str_ldf(:) = ptr_vj( -zftv(:,:,:) ) 
     251            IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( -zftv(:,:,:) ) 
     252            IF( jn == jp_sal)   str_ldf(:) = ptr_sj( -zftv(:,:,:) ) 
    253253         ENDIF 
    254254 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90

    r4990 r5682  
    2828   USE in_out_manager  ! I/O manager 
    2929   USE iom             ! I/O library 
    30 #if defined key_diaar5 
    3130   USE phycst          ! physical constants 
    3231   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    33 #endif 
    3432   USE wrk_nemo        ! Memory Allocation 
    3533   USE timing          ! Timing 
    … …  
    106104      ! 
    107105      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
     106      INTEGER  ::  ikt 
    108107      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3   ! local scalars 
    109108      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      - 
    110109      REAL(wp) ::  zcoef0, zbtr, ztra            !   -      - 
    111 #if defined key_diaar5 
    112       REAL(wp)                         ::   zztmp               ! local scalar 
    113 #endif 
    114110      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d 
    115111      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdkt, zdk1t, zdit, zdjt, ztfw  
    … …  
    149145            END DO 
    150146         END DO 
     147 
     148         ! partial cell correction 
    151149         IF( ln_zps ) THEN      ! partial steps correction at the last ocean level  
    152150            DO jj = 1, jpjm1 
    153151               DO ji = 1, fs_jpim1   ! vector opt. 
    154152! IF useless if zpshde defines pgu everywhere 
    155                   IF (mbku(ji,jj) > 1) zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
    156                   IF (mbkv(ji,jj) > 1) zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
    157                   ! (ISF) 
     153                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
     154                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
     155               END DO 
     156            END DO 
     157         ENDIF 
     158         IF( ln_zps .AND. ln_isfcav ) THEN      ! partial steps correction at the first wet level beneath a cavity 
     159            DO jj = 1, jpjm1 
     160               DO ji = 1, fs_jpim1   ! vector opt. 
    158161                  IF (miku(ji,jj) > 1) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)           
    159162                  IF (mikv(ji,jj) > 1) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)      
    160163               END DO 
    161164            END DO 
    162          ENDIF 
     165         END IF 
    163166 
    164167         !!---------------------------------------------------------------------- 
    165168         !!   II - horizontal trend  (full) 
    166169         !!---------------------------------------------------------------------- 
    167 !CDIR PARALLEL DO PRIVATE( zdk1t )  
    168          !                                                ! =============== 
    169          DO jj = 1, jpj                                 ! Horizontal slab 
    170             !                                             ! =============== 
    171             DO ji = 1, jpi   ! vector opt. 
    172                DO jk = mikt(ji,jj), jpkm1 
    173                ! 1. Vertical tracer gradient at level jk and jk+1 
    174                ! ------------------------------------------------ 
    175                ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 
    176                   zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn) - ptb(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1) 
    177                ! 
    178                   IF( jk == mikt(ji,jj) ) THEN  ;   zdkt(ji,jj,jk) = zdk1t(ji,jj,jk) 
    179                   ELSE                          ;   zdkt(ji,jj,jk) = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
    180                   ENDIF 
     170!!!!!!!!!!CDIR PARALLEL DO PRIVATE( zdk1t )  
     171            ! 1. Vertical tracer gradient at level jk and jk+1 
     172            ! ------------------------------------------------ 
     173         !  
     174         ! interior value  
     175         DO jk = 2, jpkm1                
     176            DO jj = 1, jpj 
     177               DO ji = 1, jpi   ! vector opt. 
     178                  zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn  ) - ptb(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1) 
     179                  ! 
     180                  zdkt(ji,jj,jk)  = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn  ) ) * wmask(ji,jj,jk) 
    181181               END DO 
    182182            END DO 
    183183         END DO 
    184  
    185             ! 2. Horizontal fluxes 
    186             ! --------------------    
    187          DO jj = 1 , jpjm1 
    188             DO ji = 1, fs_jpim1   ! vector opt. 
    189                DO jk = mikt(ji,jj), jpkm1 
     184         ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 
     185         zdk1t(:,:,1) = ( ptb(:,:,1,jn  ) - ptb(:,:,2,jn) ) * wmask(:,:,2) 
     186         zdkt (:,:,1) = zdk1t(:,:,1) 
     187         IF ( ln_isfcav ) THEN 
     188            DO jj = 1, jpj 
     189               DO ji = 1, jpi   ! vector opt. 
     190                  ikt = mikt(ji,jj) ! surface level 
     191                  zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn  ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1) 
     192                  zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt) 
     193               END DO 
     194            END DO 
     195         END IF 
     196 
     197         ! 2. Horizontal fluxes 
     198         ! --------------------    
     199         DO jk = 1, jpkm1 
     200            DO jj = 1 , jpjm1 
     201               DO ji = 1, fs_jpim1   ! vector opt. 
    190202                  zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) 
    191203                  zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk) 
    … …  
    208220               END DO 
    209221            END DO 
    210          END DO 
    211222 
    212223            ! II.4 Second derivative (divergence) and add to the general trend 
    213224            ! ---------------------------------------------------------------- 
    214          DO jj = 2 , jpjm1 
    215             DO ji = fs_2, fs_jpim1   ! vector opt. 
    216                DO jk = mikt(ji,jj), jpkm1 
    217                   zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     225            DO jj = 2 , jpjm1 
     226               DO ji = fs_2, fs_jpim1   ! vector opt. 
     227                  zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    218228                  ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  ) 
    219229                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
    … …  
    225235         ! 
    226236         ! "Poleward" diffusive heat or salt transports (T-S case only) 
    227          IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
     237         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
    228238            ! note sign is reversed to give down-gradient diffusive transports (#1043) 
    229             IF( jn == jp_tem)   htr_ldf(:) = ptr_vj( -zftv(:,:,:) ) 
    230             IF( jn == jp_sal)   str_ldf(:) = ptr_vj( -zftv(:,:,:) ) 
     239            IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( -zftv(:,:,:) ) 
     240            IF( jn == jp_sal)   str_ldf(:) = ptr_sj( -zftv(:,:,:) ) 
    231241         ENDIF 
    232242  
    233 #if defined key_diaar5 
    234          IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN 
    235             z2d(:,:) = 0._wp  
    236             ! note sign is reversed to give down-gradient diffusive transports (#1043) 
    237             zztmp = -1.0_wp * rau0 * rcp 
    238             DO jk = 1, jpkm1 
    239                DO jj = 2, jpjm1 
    240                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    241                      z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk)  
     243         IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN 
     244           ! 
     245           IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN 
     246               z2d(:,:) = 0._wp  
     247               DO jk = 1, jpkm1 
     248                  DO jj = 2, jpjm1 
     249                     DO ji = fs_2, fs_jpim1   ! vector opt. 
     250                        z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk)  
     251                     END DO 
    242252                  END DO 
    243253               END DO 
    244             END DO 
    245             z2d(:,:) = zztmp * z2d(:,:) 
    246             CALL lbc_lnk( z2d, 'U', -1. ) 
    247             CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction 
    248             z2d(:,:) = 0._wp  
    249             DO jk = 1, jpkm1 
    250                DO jj = 2, jpjm1 
    251                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    252                      z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk)  
     254               z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043) 
     255               CALL lbc_lnk( z2d, 'U', -1. ) 
     256               CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction 
     257               ! 
     258               z2d(:,:) = 0._wp  
     259               DO jk = 1, jpkm1 
     260                  DO jj = 2, jpjm1 
     261                     DO ji = fs_2, fs_jpim1   ! vector opt. 
     262                        z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk)  
     263                     END DO 
    253264                  END DO 
    254265               END DO 
    255             END DO 
    256             z2d(:,:) = zztmp * z2d(:,:) 
    257             CALL lbc_lnk( z2d, 'V', -1. ) 
    258             CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction 
    259          END IF 
    260 #endif 
     266               z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043) 
     267               CALL lbc_lnk( z2d, 'V', -1. ) 
     268               CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction 
     269            END IF 
     270            ! 
     271         ENDIF 
    261272 
    262273         !!---------------------------------------------------------------------- 
    … …  
    278289            DO jj = 2, jpjm1 
    279290               DO ji = fs_2, fs_jpim1   ! vector opt. 
    280                   zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
     291                  zcoef0 = - fsahtw(ji,jj,jk) * wmask(ji,jj,jk) 
    281292                  ! 
    282293                  zmsku = 1./MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)      & 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90

    r4990 r5682  
    113113      REAL(wp) ::   ze1ur, zdxt, ze2vr, ze3wr, zdyt, zdzt 
    114114      REAL(wp) ::   zah, zah_slp, zaei_slp 
    115 #if defined key_diaar5 
    116       REAL(wp) ::   zztmp              ! local scalar 
    117 #endif 
    118115      REAL(wp), POINTER, DIMENSION(:,:  ) :: z2d 
    119116      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdit, zdjt, ztfw  
    … …  
    207204      END DO 
    208205      ! 
    209 #if defined key_iomput 
    210       IF( ln_traldf_gdia .AND. cdtype == 'TRA' ) THEN 
    211          CALL wrk_alloc( jpi , jpj , jpk  , zw3d ) 
    212          DO jk=1,jpkm1 
    213             zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk)  ! u_eiv = -dpsix/dz 
    214          END DO 
    215          zw3d(:,:,jpk) = 0._wp 
    216          CALL iom_put( "uoce_eiv", zw3d )    ! i-eiv current 
    217  
    218          DO jk=1,jpk-1 
    219             zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk)  ! v_eiv = -dpsiy/dz 
    220          END DO 
    221          zw3d(:,:,jpk) = 0._wp 
    222          CALL iom_put( "voce_eiv", zw3d )    ! j-eiv current 
    223  
    224          DO jk=1,jpk-1 
    225             DO jj = 2, jpjm1 
    226                DO ji = fs_2, fs_jpim1  ! vector opt. 
    227                   zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2t(ji,jj) + & 
    228                        &    (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1t(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx 
    229                END DO 
    230             END DO 
    231          END DO 
    232          zw3d(:,:,jpk) = 0._wp 
    233          CALL iom_put( "woce_eiv", zw3d )    ! vert. eiv current 
    234          CALL wrk_dealloc( jpi , jpj , jpk  , zw3d ) 
     206      IF( iom_use("uoce_eiv") .OR. iom_use("voce_eiv") .OR. iom_use("woce_eiv") )  THEN 
     207         ! 
     208         IF( ln_traldf_gdia .AND. cdtype == 'TRA' ) THEN 
     209            CALL wrk_alloc( jpi , jpj , jpk  , zw3d ) 
     210            DO jk=1,jpkm1 
     211               zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk)  ! u_eiv = -dpsix/dz 
     212            END DO 
     213            zw3d(:,:,jpk) = 0._wp 
     214            CALL iom_put( "uoce_eiv", zw3d )    ! i-eiv current 
     215 
     216            DO jk=1,jpk-1 
     217               zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk)  ! v_eiv = -dpsiy/dz 
     218            END DO 
     219            zw3d(:,:,jpk) = 0._wp 
     220            CALL iom_put( "voce_eiv", zw3d )    ! j-eiv current 
     221 
     222            DO jk=1,jpk-1 
     223               DO jj = 2, jpjm1 
     224                  DO ji = fs_2, fs_jpim1  ! vector opt. 
     225                     zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2t(ji,jj) + & 
     226                          &    (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1t(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx 
     227                  END DO 
     228               END DO 
     229            END DO 
     230            zw3d(:,:,jpk) = 0._wp 
     231            CALL iom_put( "woce_eiv", zw3d )    ! vert. eiv current 
     232            CALL wrk_dealloc( jpi , jpj , jpk  , zw3d ) 
     233         ENDIF 
     234         ! 
    235235      ENDIF 
    236 #endif 
    237236      !                                                          ! =========== 
    238237      DO jn = 1, kjpt                                            ! tracer loop 
    … …  
    387386         ! 
    388387         !                             ! "Poleward" diffusive heat or salt transports (T-S case only) 
    389          IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
    390             IF( jn == jp_tem)   htr_ldf(:) = ptr_vj( zftv(:,:,:) )        ! 3.3  names 
    391             IF( jn == jp_sal)   str_ldf(:) = ptr_vj( zftv(:,:,:) ) 
     388         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
     389            IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( zftv(:,:,:) )        ! 3.3  names 
     390            IF( jn == jp_sal)   str_ldf(:) = ptr_sj( zftv(:,:,:) ) 
    392391         ENDIF 
    393392 
    394 #if defined key_diaar5 
    395          IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN 
    396             z2d(:,:) = 0._wp 
    397             zztmp = rau0 * rcp 
    398             DO jk = 1, jpkm1 
    399                DO jj = 2, jpjm1 
    400                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    401                      z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 
    402                   END DO 
    403                END DO 
    404             END DO 
    405             z2d(:,:) = zztmp * z2d(:,:) 
    406             CALL lbc_lnk( z2d, 'U', -1. ) 
    407             CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction 
    408             z2d(:,:) = 0._wp 
    409             DO jk = 1, jpkm1 
    410                DO jj = 2, jpjm1 
    411                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    412                      z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 
    413                   END DO 
    414                END DO 
    415             END DO 
    416             z2d(:,:) = zztmp * z2d(:,:) 
    417             CALL lbc_lnk( z2d, 'V', -1. ) 
    418             CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in j-direction 
    419          END IF 
    420 #endif 
     393         IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN 
     394           ! 
     395           IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN 
     396               z2d(:,:) = 0._wp  
     397               DO jk = 1, jpkm1 
     398                  DO jj = 2, jpjm1 
     399                     DO ji = fs_2, fs_jpim1   ! vector opt. 
     400                        z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk)  
     401                     END DO 
     402                  END DO 
     403               END DO 
     404               z2d(:,:) = rau0_rcp * z2d(:,:)  
     405               CALL lbc_lnk( z2d, 'U', -1. ) 
     406               CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction 
     407               ! 
     408               z2d(:,:) = 0._wp  
     409               DO jk = 1, jpkm1 
     410                  DO jj = 2, jpjm1 
     411                     DO ji = fs_2, fs_jpim1   ! vector opt. 
     412                        z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk)  
     413                     END DO 
     414                  END DO 
     415               END DO 
     416               z2d(:,:) = rau0_rcp * z2d(:,:)      
     417               CALL lbc_lnk( z2d, 'V', -1. ) 
     418               CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction 
     419            END IF 
     420            ! 
     421         ENDIF 
    421422         ! 
    422423      END DO 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90

    r4990 r5682  
    102102               END DO 
    103103            END DO 
    104             IF( ln_zps ) THEN      ! set gradient at partial step level 
     104            IF( ln_zps ) THEN      ! set gradient at partial step level for the last ocean cell 
    105105               DO jj = 1, jpjm1 
    106106                  DO ji = 1, fs_jpim1   ! vector opt. 
    … …  
    116116                        ztv(ji,jj,jk) = zabe2 * pgv(ji,jj,jn) 
    117117                     ENDIF 
    118                       
    119                      ! (ISH) 
     118                  END DO 
     119               END DO 
     120            ENDIF 
     121            ! (ISH) 
     122            IF( ln_zps .AND. ln_isfcav ) THEN      ! set gradient at partial step level for the first ocean cell 
     123                                                   ! into a cavity 
     124               DO jj = 1, jpjm1 
     125                  DO ji = 1, fs_jpim1   ! vector opt. 
    120126                     ! ice shelf level level MAX(2,jk) => only where ice shelf 
    121127                     iku = miku(ji,jj)  
    … …  
    148154         ! 
    149155         ! "Poleward" diffusive heat or salt transports 
    150          IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
    151             IF( jn  == jp_tem)   htr_ldf(:) = ptr_vj( ztv(:,:,:) ) 
    152             IF( jn  == jp_sal)   str_ldf(:) = ptr_vj( ztv(:,:,:) ) 
     156         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
     157            IF( jn  == jp_tem)   htr_ldf(:) = ptr_sj( ztv(:,:,:) ) 
     158            IF( jn  == jp_sal)   str_ldf(:) = ptr_sj( ztv(:,:,:) ) 
    153159         ENDIF 
    154160         !                                                  ! ================== 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/tranpc.F90

    r4990 r5682  
    99   !!            3.0  ! 2008-06  (G. Madec)  applied on ta, sa and called before tranxt in step.F90 
    1010   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA 
    11    !!            3.7  ! 2014-06  (L. Brodeau) new algorithm based on local Brunt-Vaisala freq. 
     11   !!            3.6  ! 2015-05  (L. Brodeau) new algorithm based on local Brunt-Vaisala freq. 
    1212   !!---------------------------------------------------------------------- 
    1313 
    … …  
    6464      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    6565      INTEGER  ::   inpcc        ! number of statically instable water column 
    66       INTEGER  ::   jiter, ikbot, ik, ikup, ikdown, ilayer, ikm   ! local integers 
     66      INTEGER  ::   jiter, ikbot, ikp, ikup, ikdown, ilayer, ik_low   ! local integers 
    6767      LOGICAL  ::   l_bottom_reached, l_column_treated 
    6868      REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z 
    6969      REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r2dt 
     70      REAL(wp), PARAMETER :: zn2_zero = 1.e-14_wp       ! acceptance criteria for neutrality (N2==0) 
    7071      REAL(wp), POINTER, DIMENSION(:)       ::   zvn2   ! vertical profile of N2 at 1 given point... 
    7172      REAL(wp), POINTER, DIMENSION(:,:)     ::   zvts   ! vertical profile of T and S at 1 given point... 
    … …  
    7576      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   ztrdt, ztrds   ! 3D workspace 
    7677      ! 
    77       !!LB debug: 
    78       LOGICAL, PARAMETER :: l_LB_debug = .FALSE. 
    79       INTEGER :: ilc1, jlc1, klc1, nncpu 
    80       LOGICAL :: lp_monitor_point = .FALSE. 
    81       !!LB debug. 
     78      LOGICAL, PARAMETER :: l_LB_debug = .FALSE. ! set to true if you want to follow what is 
     79      INTEGER :: ilc1, jlc1, klc1, nncpu         ! actually happening in a water column at point "ilc1, jlc1" 
     80      LOGICAL :: lp_monitor_point = .FALSE.      ! in CPU domain "nncpu" 
    8281      !!---------------------------------------------------------------------- 
    8382      ! 
    … …  
    9796         ENDIF 
    9897 
    99          !LB debug: 
    100          IF( lwp .AND. l_LB_debug ) THEN 
    101             WRITE(numout,*) 
    102             WRITE(numout,*) 'LOLO: entering tra_npc, kt, narea =', kt, narea 
    103          ENDIF 
    104          !LBdebug: Monitoring of 1 column subject to convection... 
    10598         IF( l_LB_debug ) THEN 
    106             ! Location of 1 known convection spot to follow what's happening in the water column 
    107             ilc1 = 54 ;  jlc1 = 15 ; ! Labrador ORCA1 4x4 cpus: 
    108             nncpu = 15  ; ! the CPU domain contains the convection spot 
    109             !ilc1 = 14 ;  jlc1 = 13 ; ! Labrador ORCA1 8x8 cpus: 
    110             !nncpu = 54  ; ! the CPU domain contains the convection spot 
     99            ! Location of 1 known convection site to follow what's happening in the water column 
     100            ilc1 = 45 ;  jlc1 = 3 ; !  ORCA2 4x4, Antarctic coast, more than 2 unstable portions in the  water column...            
     101            nncpu = 1  ;            ! the CPU domain contains the convection spot 
    111102            klc1 =  mbkt(ilc1,jlc1) ! bottom of the ocean for debug point...           
    112103         ENDIF 
    113          !LBdebug. 
    114  
    115          CALL eos_rab( tsa, zab )         ! after alpha and beta 
    116          CALL bn2    ( tsa, zab, zn2 )    ! after Brunt-Vaisala 
     104          
     105         CALL eos_rab( tsa, zab )         ! after alpha and beta (given on T-points) 
     106         CALL bn2    ( tsa, zab, zn2 )    ! after Brunt-Vaisala  (given on W-points) 
    117107         
    118108         inpcc = 0 
    … …  
    134124                     IF( ( ji == ilc1 ).AND.( jj == jlc1 ) ) lp_monitor_point = .TRUE. 
    135125                     ! writing only if on CPU domain where conv region is: 
    136                      lp_monitor_point = (narea == nncpu).AND.lp_monitor_point  
    137                       
    138                      IF(lp_monitor_point) THEN 
    139                         WRITE(numout,*) '' ;WRITE(numout,*) '' ; 
    140                        WRITE(numout,'("Time step = ",i6.6," !!!")')  kt 
    141                         WRITE(numout,'(" *** BEFORE anything, N^2 for point ",i3,",",i3,":" )') ji,jj 
    142                         DO jk = 1, klc1 
    143                            WRITE(numout,*) jk, zvn2(jk) 
    144                         END DO 
    145                         WRITE(numout,*) ' ' 
    146                      ENDIF 
     126                     lp_monitor_point = (narea == nncpu).AND.lp_monitor_point                       
    147127                  ENDIF                                  !LB debug  end 
    148128 
    149129                  ikbot = mbkt(ji,jj)   ! ikbot: ocean bottom T-level 
    150                   ik = 1                ! because N2 is irrelevant at the surface level (will start at ik=2) 
     130                  ikp = 1                  ! because N2 is irrelevant at the surface level (will start at ikp=2) 
    151131                  ilayer = 0 
    152132                  jiter  = 0 
    … …  
    163143                     DO WHILE ( .NOT. l_bottom_reached ) 
    164144 
    165                         ik = ik + 1 
     145                        ikp = ikp + 1 
    166146                        
    167                         !! Checking level ik for instability 
     147                        !! Testing level ikp for instability 
    168148                        !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
    169  
    170                         IF( zvn2(ik) < 0. ) THEN ! Instability found! 
    171  
    172                            ikm  = ik              ! first level whith negative N2 
    173                            ilayer = ilayer + 1    ! yet another layer found.... 
    174                            IF(jiter == 1) inpcc = inpcc + 1 
    175  
    176                            IF(l_LB_debug .AND. lp_monitor_point) & 
    177                               & WRITE(numout,*) 'Negative N2 at ik =', ikm, ' layer nb.', ilayer, & 
    178                               & ' inpcc =', inpcc 
    179  
    180                            !! Case we mix with upper regions where N2==0: 
    181                            !! All the points above ikup where N2 == 0 must also be mixed => we go 
    182                            !! upward to find a new ikup, where the layer doesn't have N2==0 
    183                            ikup = ikm 
    184                            DO jk = ikm, 2, -1 
    185                               ikup = ikup - 1 
    186                               IF( (zvn2(jk-1) > 0.).OR.(ikup == 1) ) EXIT 
    187                            END DO 
    188                            
    189                            ! adjusting ikup if the upper part of the unstable column was neutral (N2=0) 
    190                            IF((zvn2(ikup+1) == 0.).AND.(ikup /= 1)) ikup = ikup+1 ; 
    191  
    192                            
    193                            IF( lp_monitor_point )   WRITE(numout,*) ' => ikup is =', ikup, ' layer nb.', ilayer 
    194                            
     149                        IF( zvn2(ikp) <  -zn2_zero ) THEN ! Instability found! 
     150 
     151                           ilayer = ilayer + 1    ! yet another instable portion of the water column found.... 
     152 
     153                           IF( lp_monitor_point ) THEN  
     154                              WRITE(numout,*) 
     155                              IF( ilayer == 1 .AND. jiter == 1 ) THEN   ! first time a column is spoted with an instability 
     156                                 WRITE(numout,*) 
     157                                 WRITE(numout,*) 'Time step = ',kt,' !!!' 
     158                              ENDIF 
     159                              WRITE(numout,*)  ' * Iteration #',jiter,': found instable portion #',ilayer,   & 
     160                                 &                                    ' in column! Starting at ikp =', ikp 
     161                              WRITE(numout,*)  ' *** N2 for point (i,j) = ',ji,' , ',jj 
     162                              DO jk = 1, klc1 
     163                                 WRITE(numout,*) jk, zvn2(jk) 
     164                              END DO 
     165                              WRITE(numout,*) 
     166                           ENDIF 
     167                            
     168 
     169                           IF( jiter == 1 )   inpcc = inpcc + 1  
     170 
     171                           IF( lp_monitor_point )   WRITE(numout, *) 'Negative N2 at ikp =',ikp,' for layer #', ilayer 
     172 
     173                           !! ikup is the uppermost point where mixing will start: 
     174                           ikup = ikp - 1 ! ikup is always "at most at ikp-1", less if neutral levels overlying 
     175                            
     176                           !! If the points above ikp-1 have N2 == 0 they must also be mixed: 
     177                           IF( ikp > 2 ) THEN 
     178                              DO jk = ikp-1, 2, -1 
     179                                 IF( ABS(zvn2(jk)) < zn2_zero ) THEN 
     180                                    ikup = ikup - 1  ! 1 more upper level has N2=0 and must be added for the mixing 
     181                                 ELSE 
     182                                    EXIT 
     183                                 ENDIF 
     184                              END DO 
     185                           ENDIF 
     186                            
     187                           IF( ikup < 1 )   CALL ctl_stop( 'tra_npc :  PROBLEM #1') 
     188 
    195189                           zsum_temp = 0._wp 
    196190                           zsum_sali = 0._wp 
    … …  
    199193                           zsum_z    = 0._wp 
    200194                                                     
    201                            DO jk = ikup, ikbot+1      ! Inside the instable (and overlying neutral) portion of the column 
    202                               ! 
    203                               IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '     -> summing for jk =', jk 
     195                           DO jk = ikup, ikbot      ! Inside the instable (and overlying neutral) portion of the column 
    204196                              ! 
    205197                              zdz       = fse3t(ji,jj,jk) 
    … …  
    209201                              zsum_beta = zsum_beta + zvab(jk,jp_sal)*zdz 
    210202                              zsum_z    = zsum_z    + zdz 
    211                               ! 
    212                               !! EXIT if we found the bottom of the unstable portion of the water column     
    213                               IF( (zvn2(jk+1) > 0.).OR.(jk == ikbot ).OR.((jk==ikm).AND.(zvn2(jk+1) == 0.)) )   EXIT 
     203                              !                               
     204                              IF( jk == ikbot ) EXIT ! avoid array-index overshoot in case ikbot = jpk, cause we're calling jk+1 next line 
     205                              !! EXIT when we have reached the last layer that is instable (N2<0) or neutral (N2=0): 
     206                              IF( zvn2(jk+1) > zn2_zero ) EXIT 
    214207                           END DO 
    215208                           
    216                            !ik     = jk !LB remove? 
    217                            ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative N2 
    218                            
    219                            IF(l_LB_debug .AND. lp_monitor_point) & 
    220                               &    WRITE(numout,*) '  => ikdown =', ikdown, '  layer nb.', ilayer 
    221                            
    222                            ! Mixing Temperature and salinity between ikup and ikdown: 
     209                           ikdown = jk ! for the current unstable layer, ikdown is the deepest point with a negative or neutral N2 
     210                           IF( ikup == ikdown )   CALL ctl_stop( 'tra_npc :  PROBLEM #2') 
     211 
     212                           ! Mixing Temperature, salinity, alpha and beta from ikup to ikdown included: 
    223213                           zta   = zsum_temp/zsum_z 
    224214                           zsa   = zsum_sali/zsum_z 
    … …  
    226216                           zbeta = zsum_beta/zsum_z 
    227217 
    228                            IF(l_LB_debug .AND. lp_monitor_point) THEN 
     218                           IF( lp_monitor_point ) THEN 
     219                              WRITE(numout,*) 'MIXED T, S, alfa and beta between ikup =',ikup,   & 
     220                                 &            ' and ikdown =',ikdown,', in layer #',ilayer 
    229221                              WRITE(numout,*) '  => Mean temp. in that portion =', zta 
    230222                              WRITE(numout,*) '  => Mean sali. in that portion =', zsa 
    231                               WRITE(numout,*) '  => Mean Alpha in that portion =', zalfa 
     223                              WRITE(numout,*) '  => Mean Alfa in that portion =', zalfa 
    232224                              WRITE(numout,*) '  => Mean Beta  in that portion =', zbeta 
    233225                           ENDIF 
    … …  
    240232                              zvab(jk,jp_sal) = zbeta 
    241233                           END DO 
    242                            ! 
    243                            !! Before updating N2, it is possible that another unstable 
    244                            !! layer exists underneath the one we just homogeneized! 
    245                            ik = ikdown 
    246                            !  
    247                         ENDIF  ! IF( zvn2(ik+1) < 0. ) THEN 
    248                         ! 
    249                         IF( ik == ikbot ) l_bottom_reached = .TRUE. 
     234                            
     235                            
     236                           !! Updating N2 in the relvant portion of the water column 
     237                           !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion 
     238                           !! => Need to re-compute N2! will use Alpha and Beta! 
     239                            
     240                           ikup   = MAX(2,ikup)         ! ikup can never be 1 ! 
     241                           ik_low = MIN(ikdown+1,ikbot) ! we must go 1 point deeper than ikdown! 
     242                            
     243                           DO jk = ikup, ik_low              ! we must go 1 point deeper than ikdown! 
     244 
     245                              !! Interpolating alfa and beta at W point: 
     246                              zrw =  (fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk)) & 
     247                                 & / (fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk)) 
     248                              zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw 
     249                              zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw 
     250 
     251                              !! N2 at W point, doing exactly as in eosbn2.F90: 
     252                              zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     & 
     253                                 &            - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  & 
     254                                 &       / fse3w(ji,jj,jk) * tmask(ji,jj,jk) 
     255 
     256                              !! OR, faster  => just considering the vertical gradient of density 
     257                              !! as only the signa maters... 
     258                              !zvn2(jk) = (  zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     & 
     259                              !     &      - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  ) 
     260 
     261                           END DO 
     262                         
     263                           ikp = MIN(ikdown+1,ikbot) 
     264                            
     265 
     266                        ENDIF  !IF( zvn2(ikp) < 0. ) 
     267 
     268 
     269                        IF( ikp == ikbot ) l_bottom_reached = .TRUE. 
    250270                        ! 
    251271                     END DO ! DO WHILE ( .NOT. l_bottom_reached ) 
    252272 
    253                      IF( ik /= ikbot )   STOP 'ERROR: tranpc.F90 => PROBLEM #1' 
     273                     IF( ikp /= ikbot )   CALL ctl_stop( 'tra_npc :  PROBLEM #3') 
    254274                     
    255                      ! ******* At this stage ik == ikbot ! ******* 
     275                     ! ******* At this stage ikp == ikbot ! ******* 
    256276                     
    257                      IF( ilayer > 0 ) THEN 
    258                         !! least an unstable layer has been found 
    259                         !! Temperature, Salinity, Alpha and Beta have been homogenized in the unstable portion 
    260                         !! => Need to re-compute N2! will use Alpha and Beta! 
     277                     IF( ilayer > 0 ) THEN      !! least an unstable layer has been found 
    261278                        ! 
    262                         DO jk = ikup+1, ikdown+1   ! we must go 1 point deeper than ikdown!      
    263                            !! Doing exactly as in eosbn2.F90: 
    264                            !! * Except that we only are interested in the sign of N2 !!! 
    265                            !!   => just considering the vertical gradient of density 
    266                            zrw =   (fsdepw(ji,jj,jk  ) - fsdept(ji,jj,jk)) & 
    267                                & / (fsdept(ji,jj,jk-1) - fsdept(ji,jj,jk)) 
    268                            zaw = zvab(jk,jp_tem) * (1._wp - zrw) + zvab(jk-1,jp_tem) * zrw 
    269                            zbw = zvab(jk,jp_sal) * (1._wp - zrw) + zvab(jk-1,jp_sal) * zrw 
    270                            
    271                            !zvn2(jk) = grav*( zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     & 
    272                            !     &           - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )  & 
    273                            !     &       / fse3w(ji,jj,jk) * tmask(ji,jj,jk) 
    274                            zvn2(jk) = (  zaw * ( zvts(jk-1,jp_tem) - zvts(jk,jp_tem) )     & 
    275                                 &      - zbw * ( zvts(jk-1,jp_sal) - zvts(jk,jp_sal) )  )   
    276                         END DO 
    277  
    278                         IF(l_LB_debug .AND. lp_monitor_point) THEN 
    279                            WRITE(numout, '(" *** After iteration #",i3.3,", N^2 for point ",i3,",",i3,":" )') & 
    280                               & jiter, ji,jj 
     279                        IF( lp_monitor_point ) THEN 
     280                           WRITE(numout,*) 
     281                           WRITE(numout,*) 'After ',jiter,' iteration(s), we neutralized ',ilayer,' instable layer(s)' 
     282                           WRITE(numout,*) '   ==> N2 at i,j=',ji,',',jj,' now looks like this:' 
    281283                           DO jk = 1, klc1 
    282284                              WRITE(numout,*) jk, zvn2(jk) 
    283285                           END DO 
    284                            WRITE(numout,*) ' ' 
     286                           WRITE(numout,*) 
    285287                        ENDIF 
    286  
    287                         ik     = 1  ! starting again at the surface for the next iteration 
     288                        ! 
     289                        ikp    = 1     ! starting again at the surface for the next iteration 
    288290                        ilayer = 0 
    289291                     ENDIF 
    290292                     ! 
    291                      IF( ik >= ikbot ) THEN 
    292                         IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '    --- exiting jiter loop ---' 
    293                         l_column_treated = .TRUE. 
    294                      ENDIF 
     293                     IF( ikp >= ikbot )   l_column_treated = .TRUE. 
    295294                     ! 
    296295                  END DO ! DO WHILE ( .NOT. l_column_treated ) 
    … …  
    300299                  tsa(ji,jj,:,jp_sal) = zvts(:,jp_sal) 
    301300 
    302                   !! lolo:  Should we update something else???? 
    303                   !! => like alpha and beta? 
    304  
    305                   IF(l_LB_debug .AND. lp_monitor_point) WRITE(numout,*) '' 
     301                  !! LB:  Potentially some other global variable beside theta and S can be treated here 
     302                  !!      like BGC tracers. 
     303 
     304                  IF( lp_monitor_point )   WRITE(numout,*) 
    306305 
    307306               ENDIF ! IF( tmask(ji,jj,3) == 1 ) THEN 
    … …  
    321320         CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ;   CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) 
    322321         ! 
    323          IF(lwp) THEN 
    324             WRITE(numout,*) 'LOLO: exiting tra_npc, kt =', kt 
    325             WRITE(numout,*)' => number of statically instable water column : ',inpcc 
    326             WRITE(numout,*) '' ; WRITE(numout,*) '' 
     322         IF( lwp .AND. l_LB_debug ) THEN 
     323            WRITE(numout,*) 'Exiting tra_npc , kt = ',kt,', => numb. of statically instable water-columns: ', inpcc 
     324            WRITE(numout,*) 
    327325         ENDIF 
    328326         ! 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

    r4990 r5682  
    2727   USE dom_oce         ! ocean space and time domain variables  
    2828   USE sbc_oce         ! surface boundary condition: ocean 
     29   USE sbcrnf          ! river runoffs 
     30   USE sbcisf          ! ice shelf melting/freezing 
    2931   USE zdf_oce         ! ocean vertical mixing 
    3032   USE domvvl          ! variable volume 
    … …  
    4547   USE timing          ! Timing 
    4648#if defined key_agrif 
    47    USE agrif_opa_update 
    4849   USE agrif_opa_interp 
    4950#endif 
    … …  
    109110      ! Update after tracer on domain lateral boundaries 
    110111      !  
     112#if defined key_agrif 
     113      CALL Agrif_tra                     ! AGRIF zoom boundaries 
     114#endif 
     115      ! 
    111116      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1._wp )      ! local domain boundaries  (T-point, unchanged sign) 
    112117      CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1._wp ) 
    … …  
    114119#if defined key_bdy  
    115120      IF( lk_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries 
    116 #endif 
    117 #if defined key_agrif 
    118       CALL Agrif_tra                     ! AGRIF zoom boundaries 
    119121#endif 
    120122  
    … …  
    143145      ELSE                                            ! Leap-Frog + Asselin filter time stepping 
    144146         ! 
    145          IF( lk_vvl )  THEN   ;   CALL tra_nxt_vvl( kt, nit000, 'TRA', tsb, tsn, tsa, jpts )  ! variable volume level (vvl)      
    146          ELSE                 ;   CALL tra_nxt_fix( kt, nit000, 'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level  
     147         IF( lk_vvl )  THEN   ;   CALL tra_nxt_vvl( kt, nit000, rdttra, 'TRA', tsb, tsn, tsa,   & 
     148           &                                                              sbc_tsc, sbc_tsc_b, jpts )  ! variable volume level (vvl)  
     149         ELSE                 ;   CALL tra_nxt_fix( kt, nit000,         'TRA', tsb, tsn, tsa, jpts )  ! fixed    volume level  
    147150         ENDIF 
    148       ENDIF  
    149       ! 
    150 #if defined key_agrif 
    151       ! Update tracer at AGRIF zoom boundaries 
    152       IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )      ! children only 
    153 #endif       
     151      ENDIF      
    154152      ! 
    155153      ! trends computation 
    … …  
    241239 
    242240 
    243    SUBROUTINE tra_nxt_vvl( kt, kit000, cdtype, ptb, ptn, pta, kjpt ) 
     241   SUBROUTINE tra_nxt_vvl( kt, kit000, p2dt, cdtype, ptb, ptn, pta, psbc_tc, psbc_tc_b, kjpt ) 
    244242      !!---------------------------------------------------------------------- 
    245243      !!                   ***  ROUTINE tra_nxt_vvl  *** 
    … …  
    265263      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T) 
    266264      !!---------------------------------------------------------------------- 
    267       INTEGER         , INTENT(in   )                               ::   kt       ! ocean time-step index 
    268       INTEGER         , INTENT(in   )                               ::   kit000   ! first time step index 
    269       CHARACTER(len=3), INTENT(in   )                               ::   cdtype   ! =TRA or TRC (tracer indicator) 
    270       INTEGER         , INTENT(in   )                               ::   kjpt     ! number of tracers 
    271       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb      ! before tracer fields 
    272       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn      ! now tracer fields 
    273       REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta      ! tracer trend 
     265      INTEGER         , INTENT(in   )                               ::  kt       ! ocean time-step index 
     266      INTEGER         , INTENT(in   )                               ::  kit000   ! first time step index 
     267      REAL(wp)        , INTENT(in   ), DIMENSION(jpk)               ::  p2dt     ! time-step 
     268      CHARACTER(len=3), INTENT(in   )                               ::  cdtype   ! =TRA or TRC (tracer indicator) 
     269      INTEGER         , INTENT(in   )                               ::  kjpt     ! number of tracers 
     270      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptb      ! before tracer fields 
     271      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  ptn      ! now tracer fields 
     272      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::  pta      ! tracer trend 
     273      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc   ! surface tracer content 
     274      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,kjpt)      ::  psbc_tc_b ! before surface tracer content 
     275 
    274276      !!      
    275       LOGICAL  ::   ll_tra, ll_tra_hpg, ll_traqsr   ! local logical 
     277      LOGICAL  ::   ll_tra_hpg, ll_traqsr, ll_rnf, ll_isf   ! local logical 
    276278      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices 
    277279      REAL(wp) ::   zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar 
    … …  
    286288      ! 
    287289      IF( cdtype == 'TRA' )  THEN    
    288          ll_tra     = .TRUE.           ! active tracers case   
    289290         ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg 
    290291         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration 
     292         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs 
     293         IF (nn_isf .GE. 1) THEN  
     294            ll_isf = .TRUE.            ! active  tracers case  and  ice shelf melting/freezing 
     295         ELSE 
     296            ll_isf = .FALSE. 
     297         END IF 
    291298      ELSE                           
    292          ll_tra     = .FALSE.          ! passive tracers case 
    293299         ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg 
    294300         ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration 
     301         ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs 
     302         ll_isf     = .FALSE.          ! passive tracers or NO ice shelf melting/freezing 
    295303      ENDIF 
    296304      ! 
    297305      DO jn = 1, kjpt       
    298306         DO jk = 1, jpkm1 
    299             zfact1 = atfp * rdttra(jk) 
     307            zfact1 = atfp * p2dt(jk) 
    300308            zfact2 = zfact1 / rau0 
    301309            DO jj = 1, jpj 
    … …  
    315323                  ztc_f  = ztc_n  + atfp * ztc_d 
    316324                  ! 
    317                   IF( ll_tra .AND. jk == 1 ) THEN           ! first level only for T & S 
    318                       ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) ) 
    319                       ztc_f  = ztc_f  - zfact1 * ( sbc_tsc(ji,jj,jn) - sbc_tsc_b(ji,jj,jn) ) 
     325                  IF( jk == mikt(ji,jj) ) THEN           ! first level  
     326                     ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj)    - emp(ji,jj)   )  & 
     327                            &                   - (rnf_b(ji,jj)    - rnf(ji,jj)   )  & 
     328                            &                   + (fwfisf_b(ji,jj) - fwfisf(ji,jj))  ) 
     329                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) ) 
    320330                  ENDIF 
    321                   IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )   &     ! solar penetration (temperature only) 
     331 
     332                  ! solar penetration (temperature only) 
     333                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            &  
    322334                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) )  
    323335 
    324                    ze3t_f = 1.e0 / ze3t_f 
    325                    ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered 
    326                    ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta 
    327                    ! 
    328                    IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only) 
    329                       ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d ) 
    330                       pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average 
    331                    ENDIF 
     336                  ! river runoff 
     337                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          & 
     338                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) &  
     339                     &                              * fse3t_n(ji,jj,jk) / h_rnf(ji,jj) 
     340 
     341                  ! ice shelf 
     342                  IF( ll_isf ) THEN 
     343                     ! level fully include in the Losch_2008 ice shelf boundary layer 
     344                     IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          & 
     345                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  & 
     346                               &                 * fse3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) 
     347                     ! level partially include in Losch_2008 ice shelf boundary layer  
     348                     IF ( jk == misfkb(ji,jj) )                                                   & 
     349                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  & 
     350                               &                 * fse3t_n(ji,jj,jk) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj) 
     351                  END IF 
     352 
     353                  ze3t_f = 1.e0 / ze3t_f 
     354                  ptb(ji,jj,jk,jn) = ztc_f * ze3t_f       ! ptb <-- ptn filtered 
     355                  ptn(ji,jj,jk,jn) = pta(ji,jj,jk,jn)     ! ptn <-- pta 
     356                  ! 
     357                  IF( ll_tra_hpg ) THEN        ! semi-implicit hpg (T & S only) 
     358                     ze3t_d           = 1.e0   / ( ze3t_n + rbcp * ze3t_d ) 
     359                     pta(ji,jj,jk,jn) = ze3t_d * ( ztc_n  + rbcp * ztc_d  )   ! ta <-- Brown & Campana average 
     360                  ENDIF 
    332361               END DO 
    333362            END DO 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r4990 r5682  
    3232   USE wrk_nemo       ! Memory Allocation 
    3333   USE timing         ! Timing 
    34    USE sbc_ice, ONLY : lk_lim3 
    3534 
    3635   IMPLICIT NONE 
    … …  
    3837 
    3938   PUBLIC   tra_qsr       ! routine called by step.F90 (ln_traqsr=T) 
    40    PUBLIC   tra_qsr_init  ! routine called by opa.F90 
     39   PUBLIC   tra_qsr_init  ! routine called by nemogcm.F90 
    4140 
    4241   !                                 !!* Namelist namtra_qsr: penetrative solar radiation 
    … …  
    5049   REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands) 
    5150   REAL(wp), PUBLIC ::   rn_si1       !: deepest depth of extinction (water type I)       (2 bands) 
    52     
     51  
    5352   ! Module variables 
    5453   REAL(wp) ::   xsi0r                           !: inverse of rn_si0 
    … …  
    165164         CALL iom_put( 'qsr3d', etot3 )   ! Shortwave Radiation 3D distribution 
    166165         ! clem: store attenuation coefficient of the first ocean level 
    167          IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 
     166         IF ( ln_qsr_ice ) THEN 
    168167            DO jj = 1, jpj 
    169168               DO ji = 1, jpi 
    170169                  IF ( qsr(ji,jj) /= 0._wp ) THEN 
    171170                     fraqsr_1lev(ji,jj) = ( qsr_hc(ji,jj,1) / ( r1_rau0_rcp * qsr(ji,jj) ) ) 
     171                  ELSE 
     172                     fraqsr_1lev(ji,jj) = 1. 
    172173                  ENDIF 
    173174               END DO 
    … …  
    233234               END DO 
    234235               ! clem: store attenuation coefficient of the first ocean level 
    235                IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 
     236               IF ( ln_qsr_ice ) THEN 
    236237                  DO jj = 1, jpj 
    237238                     DO ji = 1, jpi 
    … …  
    256257               END DO 
    257258               ! clem: store attenuation coefficient of the first ocean level 
    258                IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 
     259               IF ( ln_qsr_ice ) THEN 
    259260                  fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp 
    260261               ENDIF 
    … …  
    279280               END DO 
    280281               ! clem: store attenuation coefficient of the first ocean level 
    281                IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 
     282               IF ( ln_qsr_ice ) THEN 
    282283                  DO jj = 1, jpj 
    283284                     DO ji = 1, jpi 
    … …  
    298299               END DO 
    299300               ! clem: store attenuation coefficient of the first ocean level 
    300                IF ( lk_lim3 .AND. ln_qsr_ice ) THEN 
     301               IF ( ln_qsr_ice ) THEN 
    301302                  fraqsr_1lev(:,:) = etot3(:,:,1) / r1_rau0_rcp 
    302303               ENDIF 
    … …  
    324325            &                    'at it= ', kt,' date= ', ndastp 
    325326         IF(lwp) WRITE(numout,*) '~~~~' 
    326          CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b', qsr_hc ) 
     327         CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b'   , qsr_hc      ) 
     328         CALL iom_rstput( kt, nitrst, numrow, 'fraqsr_1lev', fraqsr_1lev )   ! default definition in sbcssm  
    327329         ! 
    328330      ENDIF 
    … …  
    379381      ! 
    380382      IF( nn_timing == 1 )  CALL timing_start('tra_qsr_init') 
    381       ! 
    382       ! Default value for fraqsr_1lev 
    383       IF( .NOT. ln_rstart ) THEN 
    384          fraqsr_1lev(:,:) = 1._wp 
    385       ENDIF 
    386383      ! 
    387384      CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        )  
    … …  
    412409         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0 = ', rn_si0 
    413410         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1 = ', rn_si1 
    414          WRITE(numout,*) '      light penetration for ice-model LIM3     ln_qsr_ice = ', ln_qsr_ice     
    415411      ENDIF 
    416412 
    … …  
    564560      ENDIF 
    565561      ! 
     562      ! initialisation of fraqsr_1lev used in sbcssm 
     563      IF( iom_varid( numror, 'fraqsr_1lev', ldstop = .FALSE. ) > 0 ) THEN 
     564         CALL iom_get( numror, jpdom_autoglo, 'fraqsr_1lev'  , fraqsr_1lev  ) 
     565      ELSE 
     566         fraqsr_1lev(:,:) = 1._wp   ! default definition 
     567      ENDIF 
     568      ! 
    566569      CALL wrk_dealloc( jpi, jpj,      zekb, zekg, zekr        )  
    567570      CALL wrk_dealloc( jpi, jpj, jpk, ze0, ze1, ze2, ze3, zea )  

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r4990 r5682  
    99   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps 
    1010   !!             -   !  2010-09  (C. Ethe, G. Madec) Merge TRA-TRC 
     11   !!            3.6  !  2014-11  (P. Mathiot) isf melting forcing  
    1112   !!---------------------------------------------------------------------- 
    1213 
    … …  
    2021   USE sbcmod          ! ln_rnf   
    2122   USE sbcrnf          ! River runoff   
     23   USE sbcisf          ! Ice shelf    
    2224   USE traqsr          ! solar radiation penetration 
    2325   USE trd_oce         ! trends: ocean variables 
    … …  
    2628   USE in_out_manager  ! I/O manager 
    2729   USE prtctl          ! Print control 
    28    USE sbcrnf          ! River runoff   
    29    USE sbcisf          ! Ice shelf    
    30    USE sbcmod          ! ln_rnf   
    3130   USE iom 
    3231   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    … …  
    121120      REAL(wp) ::   zfact, z1_e3t, zdep 
    122121      REAL(wp) ::   zalpha, zhk 
    123       REAL(wp) ::  zt_frz, zpress 
    124122      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds 
    125123      !!---------------------------------------------------------------------- 
    … …  
    233231               DO jk = ikt, ikb - 1 
    234232               ! compute tfreez for the temperature correction (we add water at freezing temperature) 
    235 !                  zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04 
    236                   zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress ) 
    237233               ! compute trend 
    238234                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                          & 
    239                      &           + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)          & 
    240                      &               - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) & 
    241                      &           * r1_hisf_tbl(ji,jj) 
     235                     &           + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj) 
    242236                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                          & 
    243237                     &           + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) 
    … …  
    246240               ! level partially include in ice shelf boundary layer  
    247241               ! compute tfreez for the temperature correction (we add water at freezing temperature) 
    248 !               zpress = grav*rau0*fsdept(ji,jj,ikb)*1.e-04 
    249                zt_frz = -1.9 !eos_fzp( tsn(ji,jj,ikb,jp_sal), zpress ) 
    250242               ! compute trend 
    251243               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                           & 
    252                   &              + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)          & 
    253                   &                  - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) &  
    254                   &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
     244                  &              + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
    255245               tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal)                                           & 
    256246                  &              + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)  

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90

    r4990 r5682  
    8888         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    8989      END SELECT 
     90      ! DRAKKAR SSS control { 
     91      ! JMM avoid negative salinities near river outlet ! Ugly fix 
     92      ! JMM : restore negative salinities to small salinities: 
     93      WHERE ( tsa(:,:,:,jp_sal) < 0._wp )   tsa(:,:,:,jp_sal) = 0.1_wp 
    9094 
    9195      IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90

    r4990 r5682  
    122122            DO jj=1, jpj 
    123123               DO ji=1, jpi 
    124                   zwt(ji,jj,1:mikt(ji,jj)) = 0._wp 
     124                  zwt(ji,jj,1) = 0._wp 
    125125               END DO 
    126126            END DO 
    … …  
    184184            DO jj = 2, jpjm1 
    185185               DO ji = fs_2, fs_jpim1 
    186                   zwt(ji,jj,1:mikt(ji,jj)) = zwd(ji,jj,1:mikt(ji,jj)) 
    187                   DO jk = mikt(ji,jj)+1, jpkm1 
     186                  zwt(ji,jj,1) = zwd(ji,jj,1) 
     187               END DO 
     188            END DO 
     189            DO jk = 2, jpkm1 
     190               DO jj = 2, jpjm1 
     191                  DO ji = fs_2, fs_jpim1 
    188192                     zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 
    189193                  END DO 
    … …  
    196200         DO jj = 2, jpjm1 
    197201            DO ji = fs_2, fs_jpim1 
    198                ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,mikt(ji,jj)) 
    199                ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,mikt(ji,jj)) 
    200                pta(ji,jj,mikt(ji,jj),jn) = ze3tb * ptb(ji,jj,mikt(ji,jj),jn)                     & 
    201                   &                      + p2dt(mikt(ji,jj)) * ze3tn * pta(ji,jj,mikt(ji,jj),jn) 
    202                DO jk = mikt(ji,jj)+1, jpkm1 
     202               ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1) 
     203               ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,1) 
     204               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn)                     & 
     205                  &                      + p2dt(1) * ze3tn * pta(ji,jj,1,jn) 
     206            END DO 
     207         END DO 
     208         DO jk = 2, jpkm1 
     209            DO jj = 2, jpjm1 
     210               DO ji = fs_2, fs_jpim1 
    203211                  ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,jk) 
    204212                  ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t  (ji,jj,jk) 
    … …  
    213221            DO ji = fs_2, fs_jpim1 
    214222               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
    215                DO jk = jpk-2, mikt(ji,jj), -1 
     223            END DO 
     224         END DO 
     225         DO jk = jpk-2, 1, -1 
     226            DO jj = 2, jpjm1 
     227               DO ji = fs_2, fs_jpim1 
    216228                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   & 
    217229                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk) 

  • branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

    r4990 r5682  
    88   !!             -   !  2004-03  (C. Ethe)  adapted for passive tracers 
    99   !!            3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA  
     10   !!            3.6  !  2014-11  (P. Mathiot) Add zps_hde_isf (needed to open a cavity) 
    1011   !!====================================================================== 
    1112    
    … …  
    2728   PRIVATE 
    2829 
    29    PUBLIC   zps_hde    ! routine called by step.F90 
     30   PUBLIC   zps_hde     ! routine called by step.F90 
     31   PUBLIC   zps_hde_isf ! routine called by step.F90 
    3032 
    3133   !! * Substitutions 
    … …  
    4042 
    4143   SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv,   & 
     44      &                          prd, pgru, pgrv    ) 
     45      !!---------------------------------------------------------------------- 
     46      !!                     ***  ROUTINE zps_hde  *** 
     47      !!                     
     48      !! ** Purpose :   Compute the horizontal derivative of T, S and rho 
     49      !!      at u- and v-points with a linear interpolation for z-coordinate 
     50      !!      with partial steps. 
     51      !! 
     52      !! ** Method  :   In z-coord with partial steps, scale factors on last  
     53      !!      levels are different for each grid point, so that T, S and rd  
     54      !!      points are not at the same depth as in z-coord. To have horizontal 
     55      !!      gradients again, we interpolate T and S at the good depth :  
     56      !!      Linear interpolation of T, S    
     57      !!         Computation of di(tb) and dj(tb) by vertical interpolation: 
     58      !!          di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~ 
     59      !!          dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~ 
     60      !!         This formulation computes the two cases: 
     61      !!                 CASE 1                   CASE 2   
     62      !!         k-1  ___ ___________   k-1   ___ ___________ 
     63      !!                    Ti  T~                  T~  Ti+1 
     64      !!                  _____                        _____ 
     65      !!         k        |   |Ti+1     k           Ti |   | 
     66      !!                  |   |____                ____|   | 
     67      !!              ___ |   |   |           ___  |   |   | 
     68      !!                   
     69      !!      case 1->   e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then 
     70      !!          t~ = t(i+1,j  ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1) 
     71      !!        ( t~ = t(i  ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1)  ) 
     72      !!          or 
     73      !!      case 2->   e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then 
     74      !!          t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i ) 
     75      !!        ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) ) 
     76      !!          Idem for di(s) and dj(s)           
     77      !! 
     78      !!      For rho, we call eos which will compute rd~(t~,s~) at the right 
     79      !!      depth zh from interpolated T and S for the different formulations 
     80      !!      of the equation of state (eos). 
     81      !!      Gradient formulation for rho : 
     82      !!          di(rho) = rd~ - rd(i,j,k)   or   rd(i+1,j,k) - rd~ 
     83      !! 
     84      !! ** Action  : compute for top interfaces 
     85      !!              - pgtu, pgtv: horizontal gradient of tracer at u- & v-points 
     86      !!              - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points 
     87      !!---------------------------------------------------------------------- 
     88      INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index 
     89      INTEGER                              , INTENT(in   )           ::  kjpt        ! number of tracers 
     90      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields 
     91      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts  
     92      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
     93      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad of prd at u- & v-pts (bottom) 
     94      ! 
     95      INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
     96      INTEGER  ::   iku, ikv, ikum1, ikvm1   ! partial step level (ocean bottom level) at u- and v-points 
     97      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars 
     98      REAL(wp), DIMENSION(jpi,jpj)      ::  zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos 
     99      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::  zti, ztj             !  
     100      !!---------------------------------------------------------------------- 
     101      ! 
     102      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde') 
     103      ! 
     104      pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 
     105      zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 
     106      zhi (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ; 
     107      ! 
     108      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
     109         ! 
     110         DO jj = 1, jpjm1 
     111            DO ji = 1, jpim1 
     112               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
     113               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
     114               ze3wu = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
     115               ze3wv = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
     116               ! 
     117               ! i- direction 
     118               IF( ze3wu >= 0._wp ) THEN      ! case 1 
     119                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku) 
     120                  ! interpolated values of tracers 
     121                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
     122                  ! gradient of  tracers 
     123                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     124               ELSE                           ! case 2 
     125                  zmaxu = -ze3wu / fse3w(ji,jj,iku) 
     126                  ! interpolated values of tracers 
     127                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
     128                  ! gradient of tracers 
     129                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     130               ENDIF 
     131               ! 
     132               ! j- direction 
     133               IF( ze3wv >= 0._wp ) THEN      ! case 1 
     134                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv) 
     135                  ! interpolated values of tracers 
     136                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
     137                  ! gradient of tracers 
     138                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     139               ELSE                           ! case 2 
     140                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv) 
     141                  ! interpolated values of tracers 
     142                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
     143                  ! gradient of tracers 
     144                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     145               ENDIF 
     146            END DO 
     147         END DO 
     148         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     149         ! 
     150      END DO 
     151 
     152      ! horizontal derivative of density anomalies (rd) 
     153      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
     154         pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ;  
     155         DO jj = 1, jpjm1 
     156            DO ji = 1, jpim1 
     157               iku = mbku(ji,jj) 
     158               ikv = mbkv(ji,jj) 
     159               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
     160               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
     161               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)     ! i-direction: case 1 
     162               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)     ! -     -      case 2 
     163               ENDIF 
     164               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)     ! j-direction: case 1 
     165               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)     ! -     -      case 2 
     166               ENDIF 
     167            END DO 
     168         END DO 
     169 
     170         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
     171         ! step and store it in  zri, zrj for each  case 
     172         CALL eos( zti, zhi, zri )   
     173         CALL eos( ztj, zhj, zrj ) 
     174 
     175         ! Gradient of density at the last level  
     176         DO jj = 1, jpjm1 
     177            DO ji = 1, jpim1 
     178               iku = mbku(ji,jj) 
     179               ikv = mbkv(ji,jj) 
     180               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
     181               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
     182               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1 
     183               ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2 
     184               ENDIF 
     185               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1 
     186               ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj    ) )   ! j: 2 
     187               ENDIF 
     188            END DO 
     189         END DO 
     190         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions 
     191         ! 
     192      END IF 
     193      ! 
     194      IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde') 
     195      ! 
     196   END SUBROUTINE zps_hde 
     197   ! 
     198   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv,   & 
    42199      &                          prd, pgru, pgrv, pmru, pmrv, pgzu, pgzv, pge3ru, pge3rv,  & 
    43       &                   sgtu, sgtv, sgru, sgrv, smru, smrv, sgzu, sgzv, sge3ru, sge3rv ) 
     200      &                   pgtui, pgtvi, pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi ) 
    44201      !!---------------------------------------------------------------------- 
    45202      !!                     ***  ROUTINE zps_hde  *** 
    … …  
    82239      !! 
    83240      !! ** Action  : compute for top and bottom interfaces 
    84       !!              - pgtu, pgtv, sgtu, sgtv: horizontal gradient of tracer at u- & v-points 
    85       !!              - pgru, pgrv, sgru, sgtv: horizontal gradient of rho (if present) at u- & v-points 
    86       !!              - pmru, pmrv, smru, smrv: horizontal sum of rho at u- & v- point (used in dynhpg with vvl) 
    87       !!              - pgzu, pgzv, sgzu, sgzv: horizontal gradient of z at u- and v- point (used in dynhpg with vvl) 
    88       !!              - pge3ru, pge3rv, sge3ru, sge3rv: horizontal gradient of rho weighted by local e3w at u- & v-points  
     241      !!              - pgtu, pgtv, pgtui, pgtvi: horizontal gradient of tracer at u- & v-points 
     242      !!              - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points 
     243      !!              - pmru, pmrv, pmrui, pmrvi: horizontal sum of rho at u- & v- point (used in dynhpg with vvl) 
     244      !!              - pgzu, pgzv, pgzui, pgzvi: horizontal gradient of z at u- and v- point (used in dynhpg with vvl) 
     245      !!              - pge3ru, pge3rv, pge3rui, pge3rvi: horizontal gradient of rho weighted by local e3w at u- & v-points  
    89246      !!---------------------------------------------------------------------- 
    90247      INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index 
    … …  
    92249      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields 
    93250      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts  
    94       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  sgtu, sgtv  ! hor. grad. of stra at u- & v-pts (ISF) 
     251      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui, pgtvi  ! hor. grad. of stra at u- & v-pts (ISF) 
    95252      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
    96253      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv      ! hor. grad of prd at u- & v-pts (bottom) 
    … …  
    98255      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu, pgzv      ! hor. grad of z   at u- & v-pts (bottom) 
    99256      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru, pge3rv  ! hor. grad of prd weighted by local e3w at u- & v-pts (bottom) 
    100       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  sgru, sgrv      ! hor. grad of prd at u- & v-pts (top) 
    101       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  smru, smrv      ! hor. sum  of prd at u- & v-pts (top) 
    102       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  sgzu, sgzv      ! hor. grad of z   at u- & v-pts (top) 
    103       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  sge3ru, sge3rv  ! hor. grad of prd weighted by local e3w at u- & v-pts (top) 
     257      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui, pgrvi      ! hor. grad of prd at u- & v-pts (top) 
     258      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmrui, pmrvi      ! hor. sum  of prd at u- & v-pts (top) 
     259      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzui, pgzvi      ! hor. grad of z   at u- & v-pts (top) 
     260      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3rui, pge3rvi  ! hor. grad of prd weighted by local e3w at u- & v-pts (top) 
    104261      ! 
    105262      INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
    … …  
    110267      !!---------------------------------------------------------------------- 
    111268      ! 
    112       IF( nn_timing == 1 )  CALL timing_start( 'zps_hde') 
     269      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde_isf') 
    113270      ! 
    114271      pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 
    115       sgtu(:,:,:)=0.0_wp ; sgtv(:,:,:)=0.0_wp ; 
     272      pgtui(:,:,:)=0.0_wp ; pgtvi(:,:,:)=0.0_wp ; 
    116273      zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 
    117274      zhi (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ; 
    … …  
    256413                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) ) 
    257414                  ! gradient of tracers 
    258                   sgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     415                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
    259416               ELSE                           ! case 2 
    260417                  zmaxu = - ze3wu / fse3w(ji,jj,iku+1) 
    … …  
    262419                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) ) 
    263420                  ! gradient of  tracers 
    264                   sgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     421                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
    265422               ENDIF 
    266423               ! 
    … …  
    271428                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) ) 
    272429                  ! gradient of tracers 
    273                   sgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     430                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
    274431               ELSE                           ! case 2 
    275432                  zmaxv =  - ze3wv / fse3w(ji,jj,ikv+1) 
    … …  
    277434                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) ) 
    278435                  ! gradient of tracers 
    279                   sgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     436                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
    280437               ENDIF 
    281438            END DO!! 
    282439         END DO!! 
    283          CALL lbc_lnk( sgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( sgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     440         CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
    284441         ! 
    285442      END DO 
    … …  
    287444      ! horizontal derivative of density anomalies (rd) 
    288445      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
    289          sgru(:,:)  =0.0_wp ; sgrv(:,:)  =0.0_wp ; 
    290          sgzu(:,:)  =0.0_wp ; sgzv(:,:)  =0.0_wp ; 
    291          smru(:,:)  =0.0_wp ; smru(:,:)  =0.0_wp ; 
    292          sge3ru(:,:)=0.0_wp ; sge3rv(:,:)=0.0_wp ; 
     446         pgrui(:,:)  =0.0_wp ; pgrvi(:,:)  =0.0_wp ; 
     447         pgzui(:,:)  =0.0_wp ; pgzvi(:,:)  =0.0_wp ; 
     448         pmrui(:,:)  =0.0_wp ; pmrui(:,:)  =0.0_wp ; 
     449         pge3rui(:,:)=0.0_wp ; pge3rvi(:,:)=0.0_wp ; 
    293450 
    294451         DO jj = 1, jpjm1 
    … …  
    321478               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 
    322479               IF( ze3wu >= 0._wp ) THEN 
    323                  sgzu  (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku) 
    324                  sgru  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) - prd(ji,jj,iku) )          ! i: 1 
    325                  smru  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1  
    326                  sge3ru(ji,jj) = umask(ji,jj,iku+1)                                                                  & 
     480                 pgzui  (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku) 
     481                 pgrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) - prd(ji,jj,iku) )          ! i: 1 
     482                 pmrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1  
     483                 pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                  & 
    327484                                * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   & 
    328485                                   - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1 
    329486               ELSE 
    330                  sgzu  (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu) 
    331                  sgru  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) - zri(ji,jj) )      ! i: 2 
    332                  smru  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2 
    333                  sge3ru(ji,jj) = umask(ji,jj,iku+1)                                                                   & 
     487                 pgzui  (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu) 
     488                 pgrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) - zri(ji,jj) )      ! i: 2 
     489                 pmrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2 
     490                 pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                   & 
    334491                                * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  & 
    335492                                   -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2 
    336493               ENDIF 
    337494               IF( ze3wv >= 0._wp ) THEN 
    338                  sgzv  (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv)  
    339                  sgrv  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )        ! j: 1 
    340                  smrv  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1 
    341                  sge3rv(ji,jj) = vmask(ji,jj,ikv+1)                                                                  &  
     495                 pgzvi  (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv)  
     496                 pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )        ! j: 1 
     497                 pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1 
     498                 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                  &  
    342499                                * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  & 
    343500                                   - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1 
    344501                                  ! + 2 due to the formulation in density and not in anomalie in hpg sco 
    345502               ELSE 
    346                  sgzv  (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv) 
    347                  sgrv  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )     ! j: 2 
    348                  smrv  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2 
    349                  sge3rv(ji,jj) = vmask(ji,jj,ikv+1)                                                                   & 
     503                 pgzvi  (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv) 
     504                 pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )     ! j: 2 
     505                 pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2 
     506                 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                   & 
    350507                                * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) & 
    351508                                   -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2 
    … …  
    353510            END DO 
    354511         END DO 
    355          CALL lbc_lnk( sgru   , 'U', -1. )   ;   CALL lbc_lnk( sgrv   , 'V', -1. )   ! Lateral boundary conditions 
    356          CALL lbc_lnk( smru   , 'U',  1. )   ;   CALL lbc_lnk( smrv   , 'V',  1. )   ! Lateral boundary conditions 
    357          CALL lbc_lnk( sgzu   , 'U', -1. )   ;   CALL lbc_lnk( sgzv   , 'V', -1. )   ! Lateral boundary conditions 
    358          CALL lbc_lnk( sge3ru , 'U', -1. )   ;   CALL lbc_lnk( sge3rv , 'V', -1. )   ! Lateral boundary conditions 
     512         CALL lbc_lnk( pgrui   , 'U', -1. )   ;   CALL lbc_lnk( pgrvi   , 'V', -1. )   ! Lateral boundary conditions 
     513         CALL lbc_lnk( pmrui   , 'U',  1. )   ;   CALL lbc_lnk( pmrvi   , 'V',  1. )   ! Lateral boundary conditions 
     514         CALL lbc_lnk( pgzui   , 'U', -1. )   ;   CALL lbc_lnk( pgzvi   , 'V', -1. )   ! Lateral boundary conditions 
     515         CALL lbc_lnk( pge3rui , 'U', -1. )   ;   CALL lbc_lnk( pge3rvi , 'V', -1. )   ! Lateral boundary conditions 
    359516         ! 
    360517      END IF   
    361518      ! 
    362       IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde') 
    363       ! 
    364    END SUBROUTINE zps_hde 
    365  
     519      IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde_isf') 
     520      ! 
     521   END SUBROUTINE zps_hde_isf 
    366522   !!====================================================================== 
    367523END MODULE zpshde 
Note: See TracChangeset for help on using the changeset viewer.