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 10249 – NEMO

Changeset 10249


Ignore:
Timestamp:
2018-10-29T13:03:40+01:00 (6 years ago)
Author:
kingr
Message:

Merged AMM15_v3_6_STABLE_package_collate@10237

Location:
branches/UKMO/dev_r5518_AMM15_package/NEMOGCM
Files:
21 added
19 deleted
62 edited
9 copied

Legend:

Unmodified
Added
Removed

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/CONFIG/makenemo

    r10246 r10249  
    360360    cd ${NEMO_TDIR}/${NEW_CONF} || cd - 
    361361 
     362    #- 
     363    #- Get git version information if available 
     364    mkdir -p ${CONFIG_DIR}/${NEW_CONF}/BLD/inc 
     365    GITVERSIONHEADER=${CONFIG_DIR}/${NEW_CONF}/BLD/inc/gitversion.h90 
     366    echo "Extracting git commit info into" $GITVERSIONHEADER 
     367    NEMOCOMMIT=$( git describe --always --dirty ) 
     368    if [ $? -ne 0 ] 
     369    then 
     370       echo "Unable to extract git commit info!" 
     371       NEMOCOMMIT="NOT AVAILABLE" 
     372    fi 
     373    echo "#define _NEMO_COMMIT_ID_ '${NEMOCOMMIT}'" > ${GITVERSIONHEADER} 
     374    NEMOBRANCH=$( git name-rev --name-only HEAD ) 
     375    if [ $? -ne 0 ] 
     376    then 
     377       echo "Unable to extract git branch info!" 
     378       NEMOBRANCH="UNKNOWN" 
     379    fi 
     380    echo "#define _NEMO_BRANCH_ '${NEMOBRANCH}'" >> ${GITVERSIONHEADER} 
     381 
    362382#if AGRIF we do a first preprocessing 
    363383    if [ ${#x_c} -eq 0 ]; then 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90

    r10246 r10249  
    5050   USE ice 
    5151#endif 
     52   USE asminc, ONLY: ln_avgbkg 
    5253   IMPLICIT NONE 
    5354   PRIVATE 
    5455    
    5556   PUBLIC   asm_bkg_wri   !: Write out the background state 
     57 
     58  !! * variables for calculating time means 
     59   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_tavg  , sn_tavg   
     60   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   un_tavg  , vn_tavg 
     61   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   avt_tavg 
     62#if defined key_zdfgls || key_zdftke 
     63   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   en_tavg 
     64#endif 
     65   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_tavg 
     66   REAL(wp),SAVE :: numtimes_tavg     ! No of times to average over 
    5667 
    5768   !!---------------------------------------------------------------------- 
    … …  
    8192      INTEGER :: inum          ! File unit number 
    8293      REAL(wp) :: zdate        ! Date 
     94      INTEGER :: ierror 
    8395      !!----------------------------------------------------------------------- 
    8496 
    85       !                                !------------------------------------------- 
    86       IF( kt == nitbkg_r ) THEN        ! Write out background at time step nitbkg_r 
    87          !                             !-----------------------------------======== 
     97      ! If creating an averaged assim bkg, initialise on first timestep 
     98      IF ( ln_avgbkg .AND. kt == ( nn_it000 - 1) ) THEN 
     99         ! Allocate memory  
     100         ALLOCATE( tn_tavg(jpi,jpj,jpk), STAT=ierror ) 
     101         IF( ierror > 0 ) THEN 
     102            CALL ctl_stop( 'asm_wri_bkg: unable to allocate tn_tavg' )   ;   RETURN 
     103         ENDIF 
     104         tn_tavg(:,:,:)=0 
     105         ALLOCATE( sn_tavg(jpi,jpj,jpk), STAT=ierror ) 
     106         IF( ierror > 0 ) THEN 
     107            CALL ctl_stop( 'asm_wri_bkg: unable to allocate sn_tavg' )   ;   RETURN 
     108         ENDIF 
     109         sn_tavg(:,:,:)=0 
     110         ALLOCATE( un_tavg(jpi,jpj,jpk), STAT=ierror ) 
     111         IF( ierror > 0 ) THEN 
     112            CALL ctl_stop( 'asm_wri_bkg: unable to allocate un_tavg' )   ;   RETURN 
     113         ENDIF 
     114         un_tavg(:,:,:)=0 
     115         ALLOCATE( vn_tavg(jpi,jpj,jpk), STAT=ierror ) 
     116         IF( ierror > 0 ) THEN 
     117            CALL ctl_stop( 'asm_wri_bkg: unable to allocate vn_tavg' )   ;   RETURN 
     118         ENDIF 
     119         vn_tavg(:,:,:)=0 
     120         ALLOCATE( sshn_tavg(jpi,jpj), STAT=ierror ) 
     121         IF( ierror > 0 ) THEN 
     122            CALL ctl_stop( 'asm_wri_bkg: unable to allocate sshn_tavg' )   ;   RETURN 
     123         ENDIF 
     124         sshn_tavg(:,:)=0 
     125#if defined key_zdftke 
     126         ALLOCATE( en_tavg(jpi,jpj,jpk), STAT=ierror ) 
     127         IF( ierror > 0 ) THEN 
     128            CALL ctl_stop( 'asm_wri_bkg: unable to allocate en_tavg' )   ;   RETURN 
     129         ENDIF 
     130         en_tavg(:,:,:)=0 
     131#endif 
     132         ALLOCATE( avt_tavg(jpi,jpj,jpk), STAT=ierror ) 
     133         IF( ierror > 0 ) THEN 
     134            CALL ctl_stop( 'asm_wri_bkg: unable to allocate avt_tavg' )   ;   RETURN 
     135         ENDIF 
     136         avt_tavg(:,:,:)=0 
     137          
     138         numtimes_tavg = REAL ( nitavgbkg_r -  nn_it000 + 1 ) 
     139      ENDIF    
     140 
     141      ! If creating an averaged assim bkg, sum the contribution every timestep 
     142      IF ( ln_avgbkg ) THEN  
     143         IF (lwp) THEN 
     144              WRITE(numout,*) 'asm_wri_bkg : Summing assim bkg fields at timestep ',kt 
     145              WRITE(numout,*) '~~~~~~~~~~~~ ' 
     146         ENDIF 
     147 
     148         tn_tavg(:,:,:)        = tn_tavg(:,:,:) + tsn(:,:,:,jp_tem) / numtimes_tavg 
     149         sn_tavg(:,:,:)        = sn_tavg(:,:,:) + tsn(:,:,:,jp_sal) / numtimes_tavg 
     150         sshn_tavg(:,:)        = sshn_tavg(:,:) + sshn (:,:) / numtimes_tavg 
     151         un_tavg(:,:,:)        = un_tavg(:,:,:) + un(:,:,:) / numtimes_tavg 
     152         vn_tavg(:,:,:)        = vn_tavg(:,:,:) + vn(:,:,:) / numtimes_tavg 
     153         avt_tavg(:,:,:)        = avt_tavg(:,:,:) + avt(:,:,:) / numtimes_tavg 
     154#if defined key_zdftke 
     155         en_tavg(:,:,:)       = en_tavg(:,:,:) + en(:,:,:) / numtimes_tavg 
     156#endif 
     157      ENDIF 
     158      
     159 
     160      ! Write out background at time step nitbkg_r or nitavgbkg_r 
     161      IF ( ( .NOT. ln_avgbkg .AND. (kt == nitbkg_r) ) .OR. & 
     162      &          ( ln_avgbkg .AND. (kt == nitavgbkg_r) ) ) THEN 
    88163         ! 
    89164         WRITE(cl_asmbkg, FMT='(A,".nc")' ) TRIM( c_asmbkg ) 
    … …  
    97172            CALL iom_open( c_asmbkg, inum, ldwrt = .TRUE., kiolib = jprstlib) 
    98173            ! 
    99             IF( nitbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0 
    100                zdate = REAL( ndastp ) 
    101 #if defined key_zdftke 
    102                ! lk_zdftke=T :   Read turbulent kinetic energy ( en ) 
    103                IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...' 
    104                CALL tke_rst( nit000, 'READ' )               ! lk_zdftke=T :   Read turbulent kinetic energy ( en ) 
    105  
    106 #endif 
     174            ! 
     175            ! Write the information 
     176            IF ( ln_avgbkg ) THEN 
     177               IF( nitavgbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitavgbkg = 0 
     178                  zdate = REAL( ndastp ) 
     179#if defined key_zdftke 
     180                  ! lk_zdftke=T :   Read turbulent kinetic energy ( en ) 
     181                  IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...' 
     182                  CALL tke_rst( nit000, 'READ' )               ! lk_zdftke=T :   Read turbulent kinetic energy ( en ) 
     183 
     184#endif 
     185               ELSE 
     186                  zdate = REAL( ndastp ) 
     187               ENDIF 
     188               CALL iom_rstput( kt, nitavgbkg_r, inum, 'rdastp' , zdate   ) 
     189               CALL iom_rstput( kt, nitavgbkg_r, inum, 'un'     , un_tavg ) 
     190               CALL iom_rstput( kt, nitavgbkg_r, inum, 'vn'     , vn_tavg ) 
     191               CALL iom_rstput( kt, nitavgbkg_r, inum, 'tn'     , tn_tavg ) 
     192               CALL iom_rstput( kt, nitavgbkg_r, inum, 'sn'     , sn_tavg ) 
     193               CALL iom_rstput( kt, nitavgbkg_r, inum, 'sshn'   , sshn_tavg) 
     194#if defined key_zdftke 
     195               CALL iom_rstput( kt, nitavgbkg_r, inum, 'en'     , en_tavg ) 
     196#endif 
     197               CALL iom_rstput( kt, nitavgbkg_r, inum, 'avt'    , avt_tavg) 
     198               ! 
    107199            ELSE 
    108                zdate = REAL( ndastp ) 
     200               IF( nitbkg_r == nit000 - 1 ) THEN      ! Treat special case when nitbkg = 0 
     201                  zdate = REAL( ndastp ) 
     202#if defined key_zdftke 
     203                  ! lk_zdftke=T :   Read turbulent kinetic energy ( en ) 
     204                  IF(lwp) WRITE(numout,*) ' Reading TKE (en) from restart...' 
     205                  CALL tke_rst( nit000, 'READ' )               ! lk_zdftke=T :   Read turbulent kinetic energy ( en ) 
     206 
     207#endif 
     208               ELSE 
     209                  zdate = REAL( ndastp ) 
     210               ENDIF 
     211               CALL iom_rstput( kt, nitbkg_r, inum, 'rdastp' , zdate   ) 
     212               CALL iom_rstput( kt, nitbkg_r, inum, 'un'     , un                ) 
     213               CALL iom_rstput( kt, nitbkg_r, inum, 'vn'     , vn                ) 
     214               CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , tsn(:,:,:,jp_tem) ) 
     215               CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , tsn(:,:,:,jp_sal) ) 
     216               CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , sshn              ) 
     217#if defined key_zdftke 
     218               CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                ) 
     219#endif 
     220               CALL iom_rstput( kt, nitbkg_r, inum, 'avt'    , avt               ) 
     221               ! 
    109222            ENDIF 
    110             ! 
    111             !                                      ! Write the information 
    112             CALL iom_rstput( kt, nitbkg_r, inum, 'rdastp' , zdate             ) 
    113             CALL iom_rstput( kt, nitbkg_r, inum, 'un'     , un                ) 
    114             CALL iom_rstput( kt, nitbkg_r, inum, 'vn'     , vn                ) 
    115             CALL iom_rstput( kt, nitbkg_r, inum, 'tn'     , tsn(:,:,:,jp_tem) ) 
    116             CALL iom_rstput( kt, nitbkg_r, inum, 'sn'     , tsn(:,:,:,jp_sal) ) 
    117             CALL iom_rstput( kt, nitbkg_r, inum, 'sshn'   , sshn              ) 
    118 #if defined key_zdftke 
    119             CALL iom_rstput( kt, nitbkg_r, inum, 'en'     , en                ) 
    120 #endif 
    121             CALL iom_rstput( kt, nitbkg_r, inum, 'gcx'    , gcx               ) 
    122             ! 
     223             
    123224            CALL iom_close( inum ) 
     225 
    124226         ENDIF 
    125227         ! 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r10248 r10249  
    4040#endif 
    4141   USE sbc_oce          ! Surface boundary condition variables. 
     42   USE zdfmxl, ONLY :  &   
     43   &  hmld_tref,       &    
     44#if defined key_karaml 
     45   &  hmld_kara,       & 
     46   &  ln_kara,         & 
     47#endif    
     48   &  hmld,            &  
     49   &  hmlp,            & 
     50   &  hmlpt 
     51#if defined key_bdy  
     52   USE bdy_oce, ONLY: bdytmask   
     53#endif   
    4254 
    4355   IMPLICIT NONE 
    … …  
    5769#endif 
    5870   LOGICAL, PUBLIC :: ln_bkgwri = .FALSE.      !: No output of the background state fields 
     71   LOGICAL, PUBLIC :: ln_avgbkg = .FALSE.      !: No output of the mean background state fields 
    5972   LOGICAL, PUBLIC :: ln_asmiau = .FALSE.      !: No applying forcing with an assimilation increment 
    6073   LOGICAL, PUBLIC :: ln_asmdin = .FALSE.      !: No direct initialization 
    … …  
    7891   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval  
    7992   INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval 
     93   INTEGER , PUBLIC ::   nitavgbkg   !: Number of timesteps to average assim bkg [0,nitavgbkg] 
    8094   !  
    8195   INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting 
    … …  
    8599   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment 
    86100   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc 
     101 
     102   INTEGER :: mld_choice        = 4   !: choice of mld criteria to use for physics assimilation 
     103                                      !: 1) hmld      - Turbocline/mixing depth                           [W points] 
     104                                      !: 2) hmlp      - Density criterion (0.01 kg/m^3 change from 10m)   [W points] 
     105                                      !: 3) hmld_kara - Kara MLD                                          [Interpolated] 
     106                                      !: 4) hmld_tref - Temperature criterion (0.2 K change from surface) [T points] 
     107 
    87108 
    88109   !! * Substitutions 
    … …  
    117138      INTEGER :: iitiaustr_date  ! Date YYYYMMDD of IAU interval start time step 
    118139      INTEGER :: iitiaufin_date  ! Date YYYYMMDD of IAU interval final time step 
     140      INTEGER :: isurfstat       ! Local integer for status of reading surft variable 
     141      INTEGER :: iitavgbkg_date  ! Date YYYYMMDD of end of assim bkg averaging period 
    119142      ! 
    120143      REAL(wp) :: znorm        ! Normalization factor for IAU weights 
    … …  
    125148      REAL(wp) :: zdate_inc    ! Time axis in increments file 
    126149      ! 
     150      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: &  
     151          &       t_bkginc_2d  ! file for reading in 2D   
     152      !                        ! temperature increments  
     153      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: &  
     154          &       z_mld     ! Mixed layer depth  
     155           
    127156      REAL(wp), POINTER, DIMENSION(:,:) ::   hdiv   ! 2D workspace 
    128       !! 
    129       NAMELIST/nam_asminc/ ln_bkgwri,                                      & 
     157      ! 
     158      LOGICAL :: lk_surft      ! Logical: T => Increments file contains surft variable  
     159                               !               so only apply surft increments. 
     160      !! 
     161      NAMELIST/nam_asminc/ ln_bkgwri, ln_avgbkg,                           & 
    130162         &                 ln_trainc, ln_dyninc, ln_sshinc,                & 
    131163         &                 ln_asmdin, ln_asmiau,                           & 
    132164         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   & 
    133          &                 ln_salfix, salfixmin, nn_divdmp 
     165         &                 ln_salfix, salfixmin, nn_divdmp, nitavgbkg, mld_choice 
    134166      !!---------------------------------------------------------------------- 
    135167 
    … …  
    137169      ! Read Namelist nam_asminc : assimilation increment interface 
    138170      !----------------------------------------------------------------------- 
     171 
     172      ! Set default values 
     173      ln_bkgwri = .FALSE. 
     174      ln_avgbkg = .FALSE. 
     175      ln_trainc = .FALSE. 
     176      ln_dyninc = .FALSE. 
     177      ln_sshinc = .FALSE. 
    139178      ln_seaiceinc = .FALSE. 
     179      ln_asmdin = .FALSE. 
     180      ln_asmiau = .TRUE. 
     181      ln_salfix = .FALSE. 
    140182      ln_temnofreeze = .FALSE. 
     183      salfixmin = -9999 
     184      nitbkg    = 0 
     185      nitdin    = 0       
     186      nitiaustr = 1 
     187      nitiaufin = 150 
     188      niaufn    = 0 
     189      nitavgbkg = 1 
    141190 
    142191      REWIND( numnam_ref )              ! Namelist nam_asminc in reference namelist : Assimilation increment 
    … …  
    156205         WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters' 
    157206         WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri 
     207         WRITE(numout,*) '      Logical switch for writing mean background state         ln_avgbkg = ', ln_avgbkg 
    158208         WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc 
    159209         WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc 
    … …  
    166216         WRITE(numout,*) '      Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr 
    167217         WRITE(numout,*) '      Timestep of end of IAU interval in [0,nitend-nit000-1]   nitiaufin = ', nitiaufin 
     218         WRITE(numout,*) '      Number of timesteps to average assim bkg [0,nitavgbkg]   nitavgbkg = ', nitavgbkg 
    168219         WRITE(numout,*) '      Type of IAU weighting function                           niaufn    = ', niaufn 
    169220         WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix 
    170221         WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin 
     222         WRITE(numout,*) '      Choice of MLD for physics assimilation                  mld_choice = ', mld_choice 
    171223      ENDIF 
    172224 
    … …  
    175227      nitiaustr_r = nitiaustr + nit000 - 1  ! Start of IAU interval referenced to nit000 
    176228      nitiaufin_r = nitiaufin + nit000 - 1  ! End of IAU interval referenced to nit000 
     229      nitavgbkg_r = nitavgbkg + nit000 - 1  ! Averaging period referenced to nit000 
    177230 
    178231      iiauper = nitiaufin_r - nitiaustr_r + 1  ! IAU interval length 
    … …  
    184237      CALL calc_date( nit000, nitiaustr_r, ndate0, iitiaustr_date )     ! IAU start time referenced to ndate0 
    185238      CALL calc_date( nit000, nitiaufin_r, ndate0, iitiaufin_date )     ! IAU end time referenced to ndate0 
     239      CALL calc_date( nit000, nitavgbkg_r, ndate0, iitavgbkg_date )     ! End of assim bkg averaging period referenced to ndate0 
    186240      ! 
    187241      IF(lwp) THEN 
    … …  
    195249         WRITE(numout,*) '       nitiaustr_r = ', nitiaustr_r 
    196250         WRITE(numout,*) '       nitiaufin_r = ', nitiaufin_r 
     251         WRITE(numout,*) '       nitavgbkg_r = ', nitavgbkg_r 
    197252         WRITE(numout,*) 
    198253         WRITE(numout,*) '   Dates referenced to current cycle:' 
    … …  
    204259         WRITE(numout,*) '       iitiaustr_date = ', iitiaustr_date 
    205260         WRITE(numout,*) '       iitiaufin_date = ', iitiaufin_date 
     261         WRITE(numout,*) '       iitavgbkg_date = ', iitavgbkg_date 
    206262      ENDIF 
    207263 
    … …  
    246302         & CALL ctl_stop( ' nitdin :',  & 
    247303         &                ' Background time step for Direct Initialization is outside', & 
     304         &                ' the cycle interval') 
     305 
     306      IF ( nitavgbkg_r > nitend ) & 
     307         & CALL ctl_stop( ' nitavgbkg_r :',  & 
     308         &                ' Assim bkg averaging period is outside', & 
    248309         &                ' the cycle interval') 
    249310 
    … …  
    325386      !-------------------------------------------------------------------- 
    326387 
    327       ALLOCATE( t_bkginc(jpi,jpj,jpk) ) 
    328       ALLOCATE( s_bkginc(jpi,jpj,jpk) ) 
    329       ALLOCATE( u_bkginc(jpi,jpj,jpk) ) 
    330       ALLOCATE( v_bkginc(jpi,jpj,jpk) ) 
    331       ALLOCATE( ssh_bkginc(jpi,jpj)   ) 
    332       ALLOCATE( seaice_bkginc(jpi,jpj)) 
     388      IF ( ln_trainc ) THEN 
     389         ALLOCATE( t_bkginc(jpi,jpj,jpk) ) 
     390         ALLOCATE( s_bkginc(jpi,jpj,jpk) ) 
     391         t_bkginc(:,:,:) = 0.0 
     392         s_bkginc(:,:,:) = 0.0 
     393      ENDIF 
     394      IF ( ln_dyninc ) THEN  
     395         ALLOCATE( u_bkginc(jpi,jpj,jpk) ) 
     396         ALLOCATE( v_bkginc(jpi,jpj,jpk) ) 
     397         u_bkginc(:,:,:) = 0.0 
     398         v_bkginc(:,:,:) = 0.0 
     399      ENDIF 
     400      IF ( ln_sshinc ) THEN 
     401         ALLOCATE( ssh_bkginc(jpi,jpj)   ) 
     402         ssh_bkginc(:,:) = 0.0 
     403      ENDIF 
     404      IF ( ln_seaiceinc ) THEN  
     405         ALLOCATE( seaice_bkginc(jpi,jpj)) 
     406         seaice_bkginc(:,:) = 0.0 
     407      ENDIF 
    333408#if defined key_asminc 
    334409      ALLOCATE( ssh_iau(jpi,jpj)      ) 
    335 #endif 
    336       t_bkginc(:,:,:) = 0.0 
    337       s_bkginc(:,:,:) = 0.0 
    338       u_bkginc(:,:,:) = 0.0 
    339       v_bkginc(:,:,:) = 0.0 
    340       ssh_bkginc(:,:) = 0.0 
    341       seaice_bkginc(:,:) = 0.0 
    342 #if defined key_asminc 
    343410      ssh_iau(:,:)    = 0.0 
    344411#endif 
    … …  
    376443 
    377444         IF ( ln_trainc ) THEN    
    378             CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 ) 
    379             CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 ) 
    380             ! Apply the masks 
    381             t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:) 
    382             s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:) 
    383             ! Set missing increments to 0.0 rather than 1e+20 
    384             ! to allow for differences in masks 
    385             WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0 
    386             WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0 
     445 
     446            !Test if the increments file contains the surft variable. 
     447            isurfstat = iom_varid( inum, 'bckinsurft', ldstop = .FALSE. ) 
     448            IF ( isurfstat == -1 ) THEN 
     449               lk_surft = .FALSE. 
     450            ELSE 
     451               lk_surft = .TRUE. 
     452               CALL ctl_warn( ' Applying 2D temperature increment to bottom of ML: ', & 
     453            &                 ' bckinsurft found in increments file.' ) 
     454            ENDIF              
     455 
     456            IF (lk_surft) THEN  
     457                 
     458               ALLOCATE(z_mld(jpi,jpj))  
     459               SELECT CASE(mld_choice)  
     460               CASE(1)  
     461                  z_mld = hmld  
     462               CASE(2)  
     463                  z_mld = hmlp  
     464               CASE(3)  
     465#if defined key_karaml 
     466                  IF ( ln_kara ) THEN 
     467                     z_mld = hmld_kara 
     468                  ELSE 
     469                     CALL ctl_stop("Kara mixed layer not calculated as ln_kara=.false.") 
     470                  ENDIF 
     471#else 
     472                  CALL ctl_stop("Kara mixed layer not defined in current version of NEMO")  ! JW: Safety feature, should be removed 
     473                                                                                            ! once the Kara mixed layer is available  
     474#endif 
     475               CASE(4)  
     476                  z_mld = hmld_tref  
     477               END SELECT        
     478                       
     479               ALLOCATE( t_bkginc_2d(jpi,jpj) )  
     480               CALL iom_get( inum, jpdom_autoglo, 'bckinsurft', t_bkginc_2d, 1)  
     481#if defined key_bdy                 
     482               DO jk = 1,jpkm1  
     483                  WHERE( z_mld(:,:) > fsdepw(:,:,jk) )  
     484                     t_bkginc(:,:,jk) = t_bkginc_2d(:,:) * 0.5 * & 
     485                     &       ( 1 + cos( (fsdept(:,:,jk)/z_mld(:,:) ) * rpi ) ) 
     486                      
     487                     t_bkginc(:,:,jk) = t_bkginc(:,:,jk) * bdytmask(:,:)  
     488                  ELSEWHERE  
     489                     t_bkginc(:,:,jk) = 0.  
     490                  ENDWHERE  
     491               ENDDO  
     492#else  
     493               t_bkginc(:,:,:) = 0.  
     494#endif                 
     495               s_bkginc(:,:,:) = 0.  
     496                 
     497               DEALLOCATE(z_mld, t_bkginc_2d)  
     498             
     499            ELSE  
     500                
     501               CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 ) 
     502               CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 ) 
     503               ! Apply the masks 
     504               t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:) 
     505               s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:) 
     506               ! Set missing increments to 0.0 rather than 1e+20 
     507               ! to allow for differences in masks 
     508               WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0 
     509               WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0 
     510          
     511            ENDIF 
     512 
    387513         ENDIF 
    388514 
    … …  
    8901016            ENDIF 
    8911017 
     1018#if defined key_asminc 
     1019         ELSE 
     1020            ssh_iau(:,:) = 0._wp 
     1021#endif 
    8921022         ENDIF 
    8931023 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmpar.F90

    r10246 r10249  
    2929   INTEGER, PUBLIC :: nitiaufin_r   !: IAU final time step referenced to nit000 
    3030   INTEGER, PUBLIC :: nittrjfrq     !: Frequency of trajectory output for 4D-VAR 
     31   INTEGER, PUBLIC :: nitavgbkg_r   !: Averaging period for assim bkg referenced to nit000 
    3132 
    3233   !!---------------------------------------------------------------------- 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r10246 r10249  
    7171      REAL, POINTER, DIMENSION(:,:)   ::  ht_s  !: now snow thickness 
    7272#endif 
     73#if defined key_top 
     74      CHARACTER(LEN=20)                   :: cn_obc  !: type of boundary condition to apply 
     75      REAL(wp)                            :: rn_fac  !: multiplicative scaling factor 
     76      REAL(wp), POINTER, DIMENSION(:,:)   :: trc     !: now field of the tracer 
     77      LOGICAL                             :: dmp     !: obc damping term 
     78#endif 
     79 
    7380   END TYPE OBC_DATA 
    7481 
    … …  
    101108   LOGICAL, DIMENSION(jp_bdy) ::   ln_tra_dmp               !: =T Tracer damping 
    102109   LOGICAL, DIMENSION(jp_bdy) ::   ln_dyn3d_dmp             !: =T Baroclinic velocity damping 
     110    
     111!   !JT 
     112   LOGICAL, DIMENSION(jp_bdy) ::   ln_ssh_bdy               !: =T USE SSH BDY - name list switch 
     113!   !JT 
     114    
    103115   REAL(wp),    DIMENSION(jp_bdy) ::   rn_time_dmp              !: Damping time scale in days 
    104116   REAL(wp),    DIMENSION(jp_bdy) ::   rn_time_dmp_out          !: Damping time scale in days at radiation outflow points 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r10248 r10249  
    3737#endif 
    3838   USE sbcapr 
     39#if defined key_top 
     40   USE par_trc 
     41   USE trc, ONLY: trn 
     42#endif 
    3943 
    4044   IMPLICIT NONE 
    … …  
    394398#endif 
    395399      ! end jchanut tschanges 
     400       
     401       
     402      !JT use sshn (ssh now) if ln_ssh_bdy set to false in the name list 
     403      DO ib_bdy = 1, nb_bdy    
     404        nblen => idx_bdy(ib_bdy)%nblen 
     405        dta => dta_bdy(ib_bdy) 
     406          
     407        ilen1(:) = nblen(:) 
     408        !JT IF( .NOT. dta%ll_ssh ) THEN  
     409        IF( .NOT. ln_ssh_bdy(ib_bdy) ) THEN  
     410          igrd = 1 ! t Grid 
     411          DO ib = 1, ilen1(igrd) 
     412              ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     413              ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     414              dta_bdy(ib_bdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)          
     415          END DO  
     416        END IF 
     417      END DO  
    396418 
    397419      IF ( ln_apr_obc ) THEN 
    … …  
    782804            IF( dta%ll_v2d ) ALLOCATE( dta%v2d(nblen(3)) ) 
    783805         ENDIF 
    784          IF ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN 
    785             IF( dta%ll_ssh ) THEN 
    786                if(lwp) write(numout,*) '++++++ dta%ssh pointing to fnow' 
    787                jfld = jfld + 1 
    788                dta%ssh => bf(jfld)%fnow(:,1,1) 
    789             ENDIF 
     806         IF ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN          
     807            !JT  
     808            !JT allocate ssh if dta%ll_ssh set too false, as may still use it 
     809            IF (dta%ll_ssh) THEN 
     810                IF( dta%ll_ssh ) THEN 
     811                  if(lwp) write(numout,*) '++++++ dta%ssh pointing to fnow' 
     812                  jfld = jfld + 1 
     813                  dta%ssh => bf(jfld)%fnow(:,1,1) 
     814                ENDIF 
     815            ELSE 
     816              if(lwp) write(numout,*) '++++++ dta%ssh allocated space' 
     817              !ALLOCATE( dta_bdy(ib_bdy)%ssh(nblen(1)) )             
     818              ALLOCATE( dta%ssh(nblen(1)) )             
     819            ENDIF 
     820            !JT if  
     821             
    790822            IF ( dta%ll_u2d ) THEN 
    791823               IF ( ln_full_vel_array(ib_bdy) ) THEN 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90

    r10246 r10249  
    6161         CASE('zero') 
    6262            CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
     63         CASE('zerograd')  
     64            CALL bdy_dyn3d_zgrad( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )        
     65         CASE('neumann')  
     66            CALL bdy_dyn3d_nmn( idx_bdy(ib_bdy), ib_bdy ) 
    6367         CASE('orlanski') 
    6468            CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. ) 
    … …  
    117121 
    118122   END SUBROUTINE bdy_dyn3d_spe 
     123 
     124   SUBROUTINE bdy_dyn3d_zgrad( idx, dta, kt , ib_bdy )  
     125      !!----------------------------------------------------------------------  
     126      !!                  ***  SUBROUTINE bdy_dyn3d_zgrad  ***  
     127      !!  
     128      !! ** Purpose : - Enforce a zero gradient of normal velocity  
     129      !!  
     130      !!----------------------------------------------------------------------  
     131      INTEGER                     ::   kt  
     132      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices  
     133      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data  
     134      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index  
     135      !!  
     136      INTEGER  ::   jb, jk         ! dummy loop indices  
     137      INTEGER  ::   ii, ij, igrd   ! local integers  
     138      REAL(wp) ::   zwgt           ! boundary weight  
     139      INTEGER  ::   fu, fv  
     140      !!----------------------------------------------------------------------  
     141      !  
     142      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zgrad')  
     143      !  
     144      igrd = 2                      ! Copying tangential velocity into bdy points  
     145      DO jb = 1, idx%nblenrim(igrd)  
     146         DO jk = 1, jpkm1  
     147            ii   = idx%nbi(jb,igrd)  
     148            ij   = idx%nbj(jb,igrd)  
     149            fu   = ABS( ABS (NINT( idx%flagu(jb,igrd) ) ) - 1 )  
     150            ua(ii,ij,jk) = ua(ii,ij,jk) * REAL( 1 - fu) + ( ua(ii,ij+fu,jk) * umask(ii,ij+fu,jk) &  
     151                        &+ ua(ii,ij-fu,jk) * umask(ii,ij-fu,jk) ) * umask(ii,ij,jk) * REAL( fu )  
     152         END DO  
     153      END DO  
     154      !  
     155      igrd = 3                      ! Copying tangential velocity into bdy points  
     156      DO jb = 1, idx%nblenrim(igrd)  
     157         DO jk = 1, jpkm1  
     158            ii   = idx%nbi(jb,igrd)  
     159            ij   = idx%nbj(jb,igrd)  
     160            fv   = ABS( ABS (NINT( idx%flagv(jb,igrd) ) ) - 1 )  
     161            va(ii,ij,jk) = va(ii,ij,jk) * REAL( 1 - fv ) + ( va(ii+fv,ij,jk) * vmask(ii+fv,ij,jk) &  
     162                        &+ va(ii-fv,ij,jk) * vmask(ii-fv,ij,jk) ) * vmask(ii,ij,jk) * REAL( fv )  
     163         END DO  
     164      END DO  
     165      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ! Boundary points should be updated    
     166      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )     
     167      !  
     168      IF( kt .eq. nit000 ) CLOSE( unit = 102 )  
     169  
     170      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zgrad')  
     171  
     172   END SUBROUTINE bdy_dyn3d_zgrad  
    119173 
    120174   SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy ) 
    … …  
    303357   END SUBROUTINE bdy_dyn3d_dmp 
    304358 
     359   SUBROUTINE bdy_dyn3d_nmn( idx, ib_bdy )  
     360      !!----------------------------------------------------------------------  
     361      !!                 ***  SUBROUTINE bdy_dyn3d_nmn  ***  
     362      !!               
     363      !!              - Apply Neumann condition to baroclinic velocities.   
     364      !!              - Wrapper routine for bdy_nmn  
     365      !!   
     366      !!  
     367      !!----------------------------------------------------------------------  
     368      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices  
     369      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index  
     370  
     371      INTEGER  ::   jb, igrd                               ! dummy loop indices  
     372      !!----------------------------------------------------------------------  
     373  
     374      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_nmn')  
     375      !  
     376      !! Note that at this stage the ub and ua arrays contain the baroclinic velocities.   
     377      !  
     378      igrd = 2      ! Neumann bc on u-velocity;   
     379      !              
     380      CALL bdy_nmn( idx, igrd, ua )  
     381  
     382      igrd = 3      ! Neumann bc on v-velocity  
     383      !    
     384      CALL bdy_nmn( idx, igrd, va )  
     385      !  
     386      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )    ! Boundary points should be updated  
     387      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )  
     388      !  
     389      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_nmn')  
     390      !  
     391   END SUBROUTINE bdy_dyn3d_nmn  
     392  
     393 
    305394#else 
    306395   !!---------------------------------------------------------------------- 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r10248 r10249  
    105105      !! 
    106106      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend 
     107       
     108       
     109!      ! JT 
     110      NAMELIST/nambdy_ssh/ ln_ssh_bdy 
     111!      ! JT 
    107112      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    108113      !!---------------------------------------------------------------------- 
    … …  
    132137902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp ) 
    133138      IF(lwm) WRITE ( numond, nambdy ) 
     139      !JT Read nambdy_ssh namelist  
     140      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist :Unstructured open boundaries   
     141      READ  ( numnam_ref, nambdy_ssh, IOSTAT = ios, ERR = 905) 
     142905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_ssh in reference namelist', lwp ) 
     143 
     144      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries 
     145      READ  ( numnam_cfg, nambdy_ssh, IOSTAT = ios, ERR = 906) 
     146906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_ssh in configuration namelist', lwp ) 
     147      IF(lwm) WRITE ( numond, nambdy_ssh ) 
     148       
     149      IF(lwp) WRITE(numout,*) 
     150      IF(lwp) WRITE(numout,*) 'nambdy_ssh : use of ssh boundaries' 
     151      IF(lwp) WRITE(numout,*) '~~~~~~~~' 
     152      IF(lwp) WRITE(numout,*) '      ln_ssh_bdy: ' 
     153      DO ib_bdy = 1,nb_bdy 
     154        IF(lwp) WRITE(numout,*) '      ln_ssh_bdy(',ib_bdy,'): ',ln_ssh_bdy(ib_bdy) 
     155      ENDDO 
     156      IF(lwp) WRITE(numout,*) '~~~~~~~~' 
     157      IF(lwp) WRITE(numout,*)  
     158      !JT 
    134159 
    135160      ! ----------------------------------------- 
    … …  
    185210          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn2d' ) 
    186211        END SELECT 
     212         
     213        !JT override dta_bdy(ib_bdy)%ll_ssh with namelist value (ln_ssh_bdy) 
     214        IF(lwp) WRITE(numout,*) 'nambdy_ssh : use of ssh boundaries' 
     215        IF(lwp) WRITE(numout,*) '~~~~~~~~' 
     216        IF(lwp) WRITE(numout,*) '      ib_bdy: ',ib_bdy 
     217        IF(lwp) WRITE(numout,*) '      Prior to Implementation of nambdy_ssh' 
     218        IF(lwp) WRITE(numout,*) '      dta_bdy(ib_bdy)%ll_ssh: ',dta_bdy(ib_bdy)%ll_ssh 
     219         
     220        dta_bdy(ib_bdy)%ll_ssh = ln_ssh_bdy(ib_bdy) 
     221         
     222        IF(lwp) WRITE(numout,*) '      After to Implementation of nambdy_ssh' 
     223        IF(lwp) WRITE(numout,*) '      dta_bdy(ib_bdy)%ll_ssh: ',dta_bdy(ib_bdy)%ll_ssh 
     224        IF(lwp) WRITE(numout,*) '~~~~~~~~' 
     225         
     226        !JT          
     227         
    187228        IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN 
    188229           SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   !  
    … …  
    213254             dta_bdy(ib_bdy)%ll_u3d = .true. 
    214255             dta_bdy(ib_bdy)%ll_v3d = .true. 
     256          CASE('neumann')  
     257             IF(lwp) WRITE(numout,*) '      Neumann conditions'  
     258             dta_bdy(ib_bdy)%ll_u3d = .false.  
     259             dta_bdy(ib_bdy)%ll_v3d = .false.  
     260          CASE('zerograd')  
     261             IF(lwp) WRITE(numout,*) '      Zero gradient for baroclinic velocities'  
     262             dta_bdy(ib_bdy)%ll_u3d = .false.  
     263             dta_bdy(ib_bdy)%ll_v3d = .false. 
    215264          CASE('zero') 
    216265             IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
    … …  
    10871136            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    10881137               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 
    1089                idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation 
     1138               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) * 0.5 & 
     1139                                            &  *(10./ FLOAT(nn_rimwidth(ib_bdy))) ) ! JGraham:modified for rim=15 
     1140!               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation 
    10901141!               idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.  ! quadratic 
    10911142!               idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy))       ! linear 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90

    r10246 r10249  
    2626   PUBLIC   bdy_orlanski_2d     ! routine called where? 
    2727   PUBLIC   bdy_orlanski_3d     ! routine called where? 
     28   PUBLIC   bdy_nmn     ! routine called where?  
    2829 
    2930   !!---------------------------------------------------------------------- 
    … …  
    354355   END SUBROUTINE bdy_orlanski_3d 
    355356 
     357   SUBROUTINE bdy_nmn( idx, igrd, phia )  
     358      !!----------------------------------------------------------------------  
     359      !!                 ***  SUBROUTINE bdy_nmn  ***  
     360      !!                      
     361      !! ** Purpose : Duplicate the value at open boundaries, zero gradient.  
     362      !!   
     363      !!----------------------------------------------------------------------  
     364      INTEGER,                    INTENT(in)     ::   igrd     ! grid index  
     365      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   phia     ! model after 3D field (to be updated)  
     366      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices  
     367      !!   
     368      REAL(wp) ::   zcoef, zcoef1, zcoef2  
     369      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask      ! land/sea mask for field  
     370      REAL(wp), POINTER, DIMENSION(:,:)        :: bdypmask      ! land/sea mask for field  
     371      INTEGER  ::   ib, ik   ! dummy loop indices  
     372      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses  
     373      !!----------------------------------------------------------------------  
     374      !  
     375      IF( nn_timing == 1 ) CALL timing_start('bdy_nmn')  
     376      !  
     377      SELECT CASE(igrd)  
     378         CASE(1)  
     379            pmask => tmask(:,:,:)  
     380            bdypmask => bdytmask(:,:)  
     381         CASE(2)  
     382            pmask => umask(:,:,:)  
     383            bdypmask => bdyumask(:,:)  
     384         CASE(3)  
     385            pmask => vmask(:,:,:)  
     386            bdypmask => bdyvmask(:,:)  
     387         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_nmn' )  
     388      END SELECT  
     389      DO ib = 1, idx%nblenrim(igrd)  
     390         ii = idx%nbi(ib,igrd)  
     391         ij = idx%nbj(ib,igrd)  
     392         DO ik = 1, jpkm1  
     393            ! search the sense of the gradient  
     394            zcoef1 = bdypmask(ii-1,ij  )*pmask(ii-1,ij,ik) +  bdypmask(ii+1,ij  )*pmask(ii+1,ij,ik)  
     395            zcoef2 = bdypmask(ii  ,ij-1)*pmask(ii,ij-1,ik) +  bdypmask(ii  ,ij+1)*pmask(ii,ij+1,ik)  
     396            IF ( nint(zcoef1+zcoef2) == 0) THEN  
     397               ! corner **** we probably only want to set the tangentail component for the dynamics here  
     398               zcoef = pmask(ii-1,ij,ik) + pmask(ii+1,ij,ik) +  pmask(ii,ij-1,ik) +  pmask(ii,ij+1,ik)  
     399               IF (zcoef > .5_wp) THEN ! Only set none isolated points.  
     400                 phia(ii,ij,ik) = phia(ii-1,ij  ,ik) * pmask(ii-1,ij  ,ik) + &  
     401                   &              phia(ii+1,ij  ,ik) * pmask(ii+1,ij  ,ik) + &  
     402                   &              phia(ii  ,ij-1,ik) * pmask(ii  ,ij-1,ik) + &  
     403                   &              phia(ii  ,ij+1,ik) * pmask(ii  ,ij+1,ik)  
     404                 phia(ii,ij,ik) = ( phia(ii,ij,ik) / zcoef ) * pmask(ii,ij,ik)  
     405               ELSE  
     406                 phia(ii,ij,ik) = phia(ii,ij  ,ik) * pmask(ii,ij  ,ik)  
     407               ENDIF  
     408            ELSEIF ( nint(zcoef1+zcoef2) == 2) THEN  
     409               ! oblique corner **** we probably only want to set the normal component for the dynamics here  
     410               zcoef = pmask(ii-1,ij,ik)*bdypmask(ii-1,ij  ) + pmask(ii+1,ij,ik)*bdypmask(ii+1,ij  ) + &  
     411                   &   pmask(ii,ij-1,ik)*bdypmask(ii,ij -1 ) +  pmask(ii,ij+1,ik)*bdypmask(ii,ij+1  )  
     412               phia(ii,ij,ik) = phia(ii-1,ij  ,ik) * pmask(ii-1,ij  ,ik)*bdypmask(ii-1,ij  ) + &  
     413                   &            phia(ii+1,ij  ,ik) * pmask(ii+1,ij  ,ik)*bdypmask(ii+1,ij  )  + &  
     414                   &            phia(ii  ,ij-1,ik) * pmask(ii  ,ij-1,ik)*bdypmask(ii,ij -1 ) + &  
     415                   &            phia(ii  ,ij+1,ik) * pmask(ii  ,ij+1,ik)*bdypmask(ii,ij+1  )  
     416    
     417               phia(ii,ij,ik) = ( phia(ii,ij,ik) / MAX(1._wp, zcoef) ) * pmask(ii,ij,ik)  
     418            ELSE  
     419               ip = nint(bdypmask(ii+1,ij  )*pmask(ii+1,ij,ik) - bdypmask(ii-1,ij  )*pmask(ii-1,ij,ik))  
     420               jp = nint(bdypmask(ii  ,ij+1)*pmask(ii,ij+1,ik) - bdypmask(ii  ,ij-1)*pmask(ii,ij-1,ik))  
     421               phia(ii,ij,ik) = phia(ii+ip,ij+jp,ik) * pmask(ii+ip,ij+jp,ik)  
     422            ENDIF  
     423         END DO  
     424      END DO  
     425      !  
     426      IF( nn_timing == 1 ) CALL timing_stop('bdy_nmn')  
     427      !  
     428   END SUBROUTINE bdy_nmn  
     429 
    356430 
    357431#else 
    … …  
    366440      WRITE(*,*) 'bdy_orlanski_3d: You should not have seen this print! error?', kt 
    367441   END SUBROUTINE bdy_orlanski_3d 
     442   SUBROUTINE bdy_nmn( idx, igrd, phia )      ! Empty routine  
     443      WRITE(*,*) 'bdy_nmn: You should not have seen this print! error?', kt  
     444   END SUBROUTINE bdy_nmn  
    368445#endif 
    369446 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90

    r10248 r10249  
    5050      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V    (after nodal cor.) 
    5151   END TYPE TIDES_DATA 
     52   INTEGER, PUBLIC, PARAMETER                  ::   jptides_max = 15      !: Max number of tidal contituents 
     53      LOGICAL, PUBLIC                           ::   ln_harm_ana_store    !: =T Stores data for  harmonic Analysis 
     54      LOGICAL, PUBLIC                           ::   ln_harm_ana_compute     !: =T  Compute harmonic Analysis 
     55      LOGICAL, PUBLIC                           ::   ln_harmana_read         !: =T  Decide to do the analysis  
     56                                                                             !from scratch or continue previous run 
    5257 
    5358!$AGRIF_DO_NOT_TREAT 
    … …  
    9095      TYPE(MAP_POINTER), DIMENSION(jpbgrd)      ::   ibmap_ptr           !: array of pointers to nbmap 
    9196      !! 
    92       NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj 
     97      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj, ln_harm_ana_store, ln_harm_ana_compute, ln_harmana_read 
    9398      !!---------------------------------------------------------------------- 
    9499 
    … …  
    102107 
    103108      REWIND(numnam_cfg) 
     109      REWIND(numnam_ref)   ! slwa 
    104110 
    105111      DO ib_bdy = 1, nb_bdy 
    … …  
    125131            IF(lwp) WRITE(numout,*) '             assume complex conjugate   : ', ln_bdytide_conj 
    126132            IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo 
     133            IF(lwp) WRITE(numout,*) '             Use PCOMS harmonic ananalysis or not: ', ln_harm_ana_store 
     134            IF(lwp) WRITE(numout,*) '             Compute Final  harmonic ananalysis or not: ', ln_harm_ana_compute 
     135            IF(lwp) WRITE(numout,*) '             Read in previous days harmonic data or start afresh: ', ln_harmana_read 
    127136            IF(lwp) THEN  
    128137                    WRITE(numout,*) '             Tidal components: '  

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90

    r10246 r10249  
    9191      !!  
    9292      REAL(wp) ::   zwgt           ! boundary weight 
    93       INTEGER  ::   ib, ik, igrd   ! dummy loop indices 
    94       INTEGER  ::   ii, ij         ! 2D addresses 
     93      REAL(wp) ::   zcoef, zcoef1,zcoef2  
     94      INTEGER  ::   ib, ik, igrd   ! dummy loop indices  
     95      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses  
    9596      !!---------------------------------------------------------------------- 
    9697      ! 
    … …  
    160161      !!  
    161162      REAL(wp) ::   zwgt           ! boundary weight 
    162       INTEGER  ::   ib, ik, igrd   ! dummy loop indices 
    163       INTEGER  ::   ii, ij,zcoef, zcoef1,zcoef2, ip, jp   ! 2D addresses 
     163      REAL(wp) ::   zcoef, zcoef1,zcoef2  
     164      INTEGER  ::   ib, ik, igrd   ! dummy loop indices  
     165      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses  
    164166      !!---------------------------------------------------------------------- 
    165167      ! 
    … …  
    174176            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  ) 
    175177            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1) 
    176             IF ( zcoef1+zcoef2 == 0) THEN 
     178            IF ( NINT(zcoef1+zcoef2) == 0) THEN  
    177179               ! corner 
    178180               zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) +  tmask(ii,ij-1,ik) +  tmask(ii,ij+1,ik) 
    … …  
    181183                 &                    tsa(ii  ,ij-1,ik,jp_tem) * tmask(ii  ,ij-1,ik) + & 
    182184                 &                    tsa(ii  ,ij+1,ik,jp_tem) * tmask(ii  ,ij+1,ik) 
    183                tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1, zcoef) ) * tmask(ii,ij,ik) 
     185               tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1._wp, zcoef) ) * tmask(ii,ij,ik)  
    184186               tsa(ii,ij,ik,jp_sal) = tsa(ii-1,ij  ,ik,jp_sal) * tmask(ii-1,ij  ,ik) + & 
    185187                 &                    tsa(ii+1,ij  ,ik,jp_sal) * tmask(ii+1,ij  ,ik) + & 
    186188                 &                    tsa(ii  ,ij-1,ik,jp_sal) * tmask(ii  ,ij-1,ik) + & 
    187189                 &                    tsa(ii  ,ij+1,ik,jp_sal) * tmask(ii  ,ij+1,ik) 
    188                tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1, zcoef) ) * tmask(ii,ij,ik) 
     190               tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1._wp, zcoef) ) * tmask(ii,ij,ik) 
    189191            ELSE 
    190                ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  ) 
    191                jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1) 
     192               ip = NINT(bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  ))  
     193               jp = NINT(bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1))  
    192194               tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii+ip,ij+jp,ik) 
    193195               tsa(ii,ij,ik,jp_sal) = tsa(ii+ip,ij+jp,ik,jp_sal) * tmask(ii+ip,ij+jp,ik) 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90

    r10246 r10249  
    5959  PRIVATE  transport 
    6060  PRIVATE  dia_dct_wri 
     61  PRIVATE  dia_dct_wri_NOOS 
    6162 
    6263#include "domzgr_substitute.h90" 
    … …  
    6667 
    6768  !! * Module variables 
    68   INTEGER :: nn_dct        ! Frequency of computation 
    69   INTEGER :: nn_dctwri     ! Frequency of output 
    70   INTEGER :: nn_secdebug   ! Number of the section to debug 
     69  INTEGER :: nn_dct      = 1     ! Frequency of computation 
     70  INTEGER :: nn_dctwri   = 1     ! Frequency of output 
     71  INTEGER :: nn_secdebug = 0     ! Number of the section to debug 
     72  INTEGER :: nn_dct_h    = 1     ! Frequency of computation for NOOS hourly files 
     73  INTEGER :: nn_dctwri_h = 1     ! Frequency of output for NOOS hourly files 
    7174    
    72   INTEGER, PARAMETER :: nb_class_max  = 10 
    73   INTEGER, PARAMETER :: nb_sec_max    = 150 
    74   INTEGER, PARAMETER :: nb_point_max  = 2000 
    75   INTEGER, PARAMETER :: nb_type_class = 10 
    76   INTEGER, PARAMETER :: nb_3d_vars    = 3  
    77   INTEGER, PARAMETER :: nb_2d_vars    = 2  
     75  INTEGER, PARAMETER :: nb_class_max  = 12   ! maximum number of classes, i.e. depth levels or density classes 
     76  INTEGER, PARAMETER :: nb_sec_max    = 30   ! maximum number of sections 
     77  INTEGER, PARAMETER :: nb_point_max  = 300  ! maximum number of points in a single section 
     78  INTEGER, PARAMETER :: nb_type_class       = 14   ! types of calculations, i.e. pos transport, neg transport, heat transport, salt transport 
     79  INTEGER, PARAMETER :: nb_3d_vars    = 5 
     80  INTEGER, PARAMETER :: nb_2d_vars    = 2 
    7881  INTEGER            :: nb_sec  
    7982 
    … …  
    102105                                                     zlay              ! level classes      (99 if you don't want) 
    103106     REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport     ! transport output 
     107     REAL(wp), DIMENSION(nb_type_class,nb_class_max)  :: transport_h   ! transport output 
    104108     REAL(wp)                                         :: slopeSection  ! slope of the section 
    105109     INTEGER                                          :: nb_point      ! number of points in the section 
    … …  
    111115  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d  
    112116  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d   
     117  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d_h 
     118  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d_h 
     119  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  z_hr_output 
    113120 
    114121   !! $Id$ 
    … …  
    120127     !!                   ***  FUNCTION diadct_alloc  ***  
    121128     !!----------------------------------------------------------------------  
    122      INTEGER :: ierr(2)  
     129     INTEGER :: ierr(5)  
    123130     !!----------------------------------------------------------------------  
    124131 
    125      ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(1) )  
    126      ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=ierr(2) )  
     132     ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk)  , STAT=ierr(1) )  
     133     ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)      , STAT=ierr(2) )  
     134     ALLOCATE(transports_3d_h(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(3) ) 
     135     ALLOCATE(transports_2d_h(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=ierr(4) ) 
     136     ALLOCATE(z_hr_output(nb_sec_max,168,nb_class_max)               , STAT=ierr(5) ) ! 168 = 24 * 7days 
    127137 
    128138     diadct_alloc = MAXVAL( ierr )  
    … …  
    153163     IF(lwm) WRITE ( numond, namdct ) 
    154164 
     165     IF( ln_NOOS ) THEN 
     166        nn_dct=3600./rdt         ! hard coded for NOOS transects, to give 25 hour means  
     167        nn_dctwri=86400./rdt 
     168        nn_dct_h=1       ! hard coded for NOOS transects, to give hourly data 
     169        nn_dctwri_h=3600./rdt 
     170     ENDIF 
     171 
    155172     IF( lwp ) THEN 
    156173        WRITE(numout,*) " " 
    157174        WRITE(numout,*) "diadct_init: compute transports through sections " 
    158175        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~" 
    159         WRITE(numout,*) "       Frequency of computation: nn_dct    = ",nn_dct 
    160         WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri 
     176        IF( ln_NOOS ) THEN 
     177           WRITE(numout,*) "       Frequency of computation hard coded to be every hour: nn_dct    = ",nn_dct 
     178           WRITE(numout,*) "       Frequency of write hard coded to average 25 instantaneous hour values: nn_dctwri = ",nn_dctwri 
     179           WRITE(numout,*) "       Frequency of hourly computation hard coded to be every timestep: nn_dct_h  = ",nn_dct_h 
     180           WRITE(numout,*) "       Frequency of hourly write hard coded to every hour: nn_dctwri_h = ",nn_dctwri_h 
     181        ELSE 
     182           WRITE(numout,*) "       Frequency of computation: nn_dct    = ",nn_dct 
     183           WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri 
     184        ENDIF 
    161185 
    162186        IF      ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN 
    … …  
    177201     !open output file 
    178202     IF( lwm ) THEN 
    179         CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
    180         CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
    181         CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     203        IF( ln_NOOS ) THEN 
     204           CALL ctl_opn( numdct_NOOS  ,'NOOS_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     205           CALL ctl_opn( numdct_NOOS_h,'NOOS_transport_h', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     206        ELSE 
     207           CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     208           CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     209           CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     210        ENDIF 
    182211     ENDIF 
    183212 
    184213     ! Initialise arrays to zero  
    185      transports_3d(:,:,:,:)=0.0  
    186      transports_2d(:,:,:)  =0.0  
     214     transports_3d(:,:,:,:)  =0.0  
     215     transports_2d(:,:,:)    =0.0  
     216     transports_3d_h(:,:,:,:)=0._wp 
     217     transports_2d_h(:,:,:)  =0._wp 
     218     z_hr_output(:,:,:)      =0._wp 
    187219 
    188220     IF( nn_timing == 1 )   CALL timing_stop('dia_dct_init') 
    … …  
    229261        CALL wrk_alloc( nb_sec_max,nb_type_class,nb_class_max , zsum  ) 
    230262     ENDIF     
    231   
     263     lldebug=.TRUE. 
    232264     ! Initialise arrays 
    233265     zwork(:) = 0.0  
    … …  
    243275  
    244276     ! Compute transport and write only at nn_dctwri 
    245      IF( MOD(kt,nn_dct)==0 ) THEN  
     277     IF( MOD(kt,nn_dct)==0 .or. &                ! compute transport every nn_dct time steps 
     278         (ln_NOOS .and. kt==nn_it000 ) )  THEN   ! also include first time step when calculating NOOS 25 hour averages 
    246279 
    247280        DO jsec=1,nb_sec 
    248281 
    249282           !debug this section computing ? 
    250            lldebug=.FALSE. 
    251283           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct-1 .AND. lwp ) lldebug=.TRUE.  
    252284 
    … …  
    271303           !Sum on all procs  
    272304           IF( lk_mpp )THEN 
     305              zsum(:,:,:)=0.0_wp 
    273306              ish(1)  =  nb_sec_max*nb_type_class*nb_class_max  
    274307              ish2    = (/nb_sec_max,nb_type_class,nb_class_max/) 
    … …  
    283316           DO jsec=1,nb_sec 
    284317 
    285               IF( lwm )CALL dia_dct_wri(kt,jsec,secs(jsec)) 
     318              IF( lwm .and. .not. ln_NOOS )CALL dia_dct_wri(kt,jsec,secs(jsec)) 
     319              IF( lwm .and.       ln_NOOS )CALL dia_dct_wri_NOOS(kt,jsec,secs(jsec))   ! use NOOS specific formatting 
    286320             
    287321              !nullify transports values after writing 
    288322              transports_3d(:,jsec,:,:)=0. 
    289323              transports_2d(:,jsec,:  )=0. 
    290               secs(jsec)%transport(:,:)=0.   
     324              secs(jsec)%transport(:,:)=0.  
     325              IF ( ln_NOOS ) CALL transport(secs(jsec),lldebug,jsec)  ! reinitialise for next 25 hour instantaneous average (overlapping values) 
    291326 
    292327           ENDDO 
    … …  
    295330 
    296331     ENDIF 
     332 
     333     IF ( MOD(kt,nn_dct_h)==0 ) THEN            ! compute transport every nn_dct_h time steps 
     334 
     335        DO jsec=1,nb_sec 
     336 
     337           !lldebug=.FALSE. 
     338           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct_h-1 .AND. lwp ) lldebug=.TRUE.  
     339 
     340           !Compute transport through section   
     341           CALL transport_h(secs(jsec),lldebug,jsec)  
     342 
     343        ENDDO 
     344              
     345        IF( MOD(kt,nn_dctwri_h)==0 )THEN 
     346 
     347           IF( lwp .AND. kt==nit000+nn_dctwri_h-1 )WRITE(numout,*)"      diadct: average and write hourly files at kt = ",kt          
     348   
     349           !! divide arrays by nn_dctwri/nn_dct to obtain average 
     350           transports_3d_h(:,:,:,:)=transports_3d_h(:,:,:,:)/(nn_dctwri_h/nn_dct_h) 
     351           transports_2d_h(:,:,:)  =transports_2d_h(:,:,:)  /(nn_dctwri_h/nn_dct_h) 
     352 
     353           ! Sum over each class 
     354           DO jsec=1,nb_sec 
     355              CALL dia_dct_sum_h(secs(jsec),jsec) 
     356           ENDDO 
     357  
     358           !Sum on all procs  
     359          IF( lk_mpp )THEN 
     360              ish(1)  =  nb_sec_max*nb_type_class*nb_class_max  
     361              ish2    = (/nb_sec_max,nb_type_class,nb_class_max/) 
     362              DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport_h(:,:) ; ENDDO 
     363              zwork(:)= RESHAPE(zsum(:,:,:), ish ) 
     364              CALL mpp_sum(zwork, ish(1)) 
     365              zsum(:,:,:)= RESHAPE(zwork,ish2) 
     366              DO jsec=1,nb_sec ; secs(jsec)%transport_h(:,:) = zsum(jsec,:,:) ; ENDDO 
     367           ENDIF 
     368 
     369           !Write the transport 
     370           DO jsec=1,nb_sec 
     371 
     372              IF( lwp .and.       ln_NOOS )CALL dia_dct_wri_NOOS_h(kt/nn_dctwri_h,jsec,secs(jsec))   ! use NOOS specific formatting 
     373             
     374              !nullify transports values after writing 
     375              transports_3d_h(:,jsec,:,:)=0.0 
     376              transports_2d_h(:,jsec,:)=0.0 
     377              secs(jsec)%transport_h(:,:)=0.   
     378              IF ( ln_NOOS ) CALL transport_h(secs(jsec),lldebug,jsec)  ! reinitialise for next 25 hour instantaneous average (overlapping values) 
     379 
     380           ENDDO 
     381 
     382        ENDIF  
     383 
     384     ENDIF     
    297385 
    298386     IF( lk_mpp )THEN 
    … …  
    356444        secs(jsec)%zlay         = 99._wp          
    357445        secs(jsec)%transport    =  0._wp   ; secs(jsec)%nb_point       = 0 
     446        secs(jsec)%transport_h  =  0._wp   ; secs(jsec)%nb_point       = 0 
    358447 
    359448        !read section's number / name / computing choices / classes / slopeSection / points number 
    360449        !----------------------------------------------------------------------------------------- 
    361450        READ(numdct_in,iostat=iost)isec 
    362         IF (iost .NE. 0 )EXIT !end of file  
     451        IF (iost .NE. 0 )EXIT !end of file 
    363452        WRITE(cltmp,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' and jsec= ',jsec 
    364         IF( jsec .NE. isec )  CALL ctl_stop( cltmp ) 
     453        IF( jsec .NE. isec )CALL ctl_stop( cltmp ) 
    365454 
    366455        IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )WRITE(numout,*)"isec ",isec  
    … …  
    380469        READ(numdct_in)iptglo 
    381470 
     471        IF ( ln_NOOS ) THEN 
     472           WRITE(numout,*) 'Section name = ',TRIM(secs(jsec)%name) 
     473           WRITE(numout,*) 'coordSec = ',secs(jsec)%coordSec 
     474           WRITE(numout,*) 'iptglo = ',iptglo 
     475        ENDIF 
     476 
    382477        !debug 
    383478        !----- 
    … …  
    405500           !read points'coordinates and directions  
    406501           !-------------------------------------- 
     502           IF ( ln_NOOS ) THEN 
     503              WRITE(numout,*) 'Coords and direction read' 
     504           ENDIF 
     505 
    407506           coordtemp(:) = POINT_SECTION(0,0) !list of points read 
    408507           directemp(:) = 0                  !value of directions of each points 
    … …  
    422521              ENDDO                   
    423522           ENDIF 
    424   
     523 
    425524           !Now each proc selects only points that are in its domain: 
    426525           !-------------------------------------------------------- 
    … …  
    624723     !!-------------------------------------------------------- 
    625724 
    626      IF( ld_debug )WRITE(numout,*)'      Compute transport' 
     725     !!NIALL IF( ld_debug )WRITE(numout,*)'      Compute transport' 
    627726 
    628727     !---------------------------! 
    … …  
    710809              SELECT CASE( sec%direction(jseg) )  
    711810              CASE(0,1)  
    712                  ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )  
    713                  zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )  
    714                  zrhop = interp(k%I,k%J,jk,'V',rhop)  
    715                  zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)  
     811                 ztn   = interp(k%I,k%J,jk,'V',0 )  
     812                 zsn   = interp(k%I,k%J,jk,'V',1 )  
     813                 zrhop = interp(k%I,k%J,jk,'V',2 )  
     814                 zrhoi = interp(k%I,k%J,jk,'V',3 )  
    716815                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)  
    717816              CASE(2,3)  
    718                  ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )  
    719                  zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )  
    720                  zrhop = interp(k%I,k%J,jk,'U',rhop)  
    721                  zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)  
     817                 ztn   = interp(k%I,k%J,jk,'U',0 )  
     818                 zsn   = interp(k%I,k%J,jk,'U',1 )  
     819                 zrhop = interp(k%I,k%J,jk,'U',2 )  
     820                 zrhoi = interp(k%I,k%J,jk,'U',3 )  
    722821                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)   
    723822              END SELECT  
    … …  
    754853                 transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk)  + zTnorm * ztn * zrhop * rcp 
    755854                 transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk)  + zTnorm * zsn * zrhop * 0.001 
     855                 transports_3d(4,jsec,jseg,jk) = transports_3d(4,jsec,jseg,jk)  + zTnorm * 4.e+3_wp * (ztn+273.15) * 1026._wp 
     856                 transports_3d(5,jsec,jseg,jk) = transports_3d(5,jsec,jseg,jk)  + zTnorm * 0.001 * zsn * 1026._wp 
    756857              ENDIF 
    757858    
    … …  
    809910  END SUBROUTINE transport 
    810911   
     912  SUBROUTINE transport_h(sec,ld_debug,jsec) 
     913     !!------------------------------------------------------------------------------------------- 
     914     !!                     ***  ROUTINE hourly_transport  *** 
     915     !! 
     916     !!              Exactly as routine transport but for data summed at 
     917     !!              each time step and averaged each hour 
     918     !! 
     919     !!  Purpose ::  Compute the transport for each point in a section 
     920     !! 
     921     !!  Method  ::  Loop over each segment, and each vertical level and add the transport 
     922     !!              Be aware :           
     923     !!              One section is a sum of segments 
     924     !!              One segment is defined by 2 consecutive points in sec%listPoint 
     925     !!              All points of sec%listPoint are positioned on the F-point of the cell 
     926     !!  
     927     !!              There are two loops:                  
     928     !!              loop on the segment between 2 nodes 
     929     !!              loop on the level jk 
     930     !! 
     931     !! ** Output: Arrays containing the volume,density,salinity,temperature etc 
     932     !!            transports for each point in a section, summed over each nn_dct. 
     933     !! 
     934     !!------------------------------------------------------------------------------------------- 
     935     !! * Arguments 
     936     TYPE(SECTION),INTENT(INOUT) :: sec 
     937     LOGICAL      ,INTENT(IN)    :: ld_debug 
     938     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section 
     939     
     940     !! * Local variables 
     941     INTEGER             :: jk,jseg,jclass,    &!loop on level/segment/classes  
     942                            isgnu  , isgnv      ! 
     943     REAL(wp):: zumid        , zvmid        ,  &!U/V velocity on a cell segment 
     944                zumid_ice    , zvmid_ice    ,  &!U/V ice velocity 
     945                zTnorm                          !transport of velocity through one cell's sides 
     946     REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 
     947 
     948     TYPE(POINT_SECTION) :: k 
     949     !!-------------------------------------------------------- 
     950 
     951     !!NIALL IF( ld_debug )WRITE(numout,*)'      Compute transport' 
     952 
     953     !---------------------------! 
     954     !  COMPUTE TRANSPORT        ! 
     955     !---------------------------! 
     956     IF(sec%nb_point .NE. 0)THEN    
     957 
     958        !---------------------------------------------------------------------------------------------------- 
     959        !Compute sign for velocities: 
     960        ! 
     961        !convention: 
     962        !   non horizontal section: direction + is toward left hand of section 
     963        !       horizontal section: direction + is toward north of section 
     964        ! 
     965        ! 
     966        !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0 
     967        !       ----------------      -----------------     ---------------             -------------- 
     968        ! 
     969        !   isgnv=1         direction +       
     970        !  ______         _____             ______                                                    
     971        !        |           //|            |                  |                         direction +    
     972        !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\ 
     973        !        |_______  //         ______|    \\            | ---\                        | 
     974        !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________ 
     975        !               |             |          __\\|         |                     
     976        !               |             |     direction +        |                      isgnv=1                                  
     977        !                                                       
     978        !---------------------------------------------------------------------------------------------------- 
     979        isgnu = 1 
     980        IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1  
     981        ELSE                                ; isgnv =  1 
     982        ENDIF 
     983 
     984        IF( ld_debug )write(numout,*)"isgnu isgnv ",isgnu,isgnv 
     985 
     986        !--------------------------------------! 
     987        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  ! 
     988        !--------------------------------------! 
     989        DO jseg=1,MAX(sec%nb_point-1,0) 
     990               
     991           !------------------------------------------------------------------------------------------- 
     992           ! Select the appropriate coordinate for computing the velocity of the segment 
     993           ! 
     994           !                      CASE(0)                                    Case (2) 
     995           !                      -------                                    -------- 
     996           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)       
     997           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               | 
     998           !                                                                            | 
     999           !                                                                            | 
     1000           !                                                                            | 
     1001           !                      Case (3)                                            U(i,j) 
     1002           !                      --------                                              | 
     1003           !                                                                            | 
     1004           !  listPoint(jseg+1) F(i,j+1)                                                | 
     1005           !                        |                                                   | 
     1006           !                        |                                                   | 
     1007           !                        |                                 listPoint(jseg+1) F(i,j-1) 
     1008           !                        |                                             
     1009           !                        |                                             
     1010           !                     U(i,j+1)                                             
     1011           !                        |                                       Case(1)      
     1012           !                        |                                       ------       
     1013           !                        |                                             
     1014           !                        |                 listPoint(jseg+1)             listPoint(jseg)                            
     1015           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                            
     1016           ! listPoint(jseg)     F(i,j) 
     1017           !  
     1018           !------------------------------------------------------------------------------------------- 
     1019 
     1020           SELECT CASE( sec%direction(jseg) ) 
     1021           CASE(0)  ;   k = sec%listPoint(jseg) 
     1022           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 
     1023           CASE(2)  ;   k = sec%listPoint(jseg) 
     1024           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 
     1025           END SELECT 
     1026 
     1027           !---------------------------| 
     1028           !     LOOP ON THE LEVEL     | 
     1029           !---------------------------| 
     1030           !Sum of the transport on the vertical  
     1031           DO jk=1,mbathy(k%I,k%J) 
     1032 
     1033              ! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point 
     1034              SELECT CASE( sec%direction(jseg) ) 
     1035              CASE(0,1) 
     1036                 ztn   = interp(k%I,k%J,jk,'V',0 ) 
     1037                 zsn   = interp(k%I,k%J,jk,'V',1) 
     1038                 zrhop = interp(k%I,k%J,jk,'V',2) 
     1039                 zrhoi = interp(k%I,k%J,jk,'V',3) 
     1040                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
     1041              CASE(2,3) 
     1042                 ztn   = interp(k%I,k%J,jk,'U',0) 
     1043                 zsn   = interp(k%I,k%J,jk,'U',1) 
     1044                 zrhop = interp(k%I,k%J,jk,'U',2) 
     1045                 zrhoi = interp(k%I,k%J,jk,'U',3) 
     1046                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
     1047              END SELECT 
     1048 
     1049              zfsdep= fsdept(k%I,k%J,jk) 
     1050  
     1051              !compute velocity with the correct direction 
     1052              SELECT CASE( sec%direction(jseg) ) 
     1053              CASE(0,1)   
     1054                 zumid=0. 
     1055                 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 
     1056              CASE(2,3) 
     1057                 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 
     1058                 zvmid=0. 
     1059              END SELECT 
     1060 
     1061              !zTnorm=transport through one cell; 
     1062              !velocity* cell's length * cell's thickness 
     1063              zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     & 
     1064                     zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk) 
     1065 
     1066#if ! defined key_vvl 
     1067              !add transport due to free surface 
     1068              IF( jk==1 )THEN 
     1069                 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + & 
     1070                                   zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 
     1071              ENDIF 
     1072#endif 
     1073              !COMPUTE TRANSPORT  
     1074 
     1075              transports_3d_h(1,jsec,jseg,jk) = transports_3d_h(1,jsec,jseg,jk) + zTnorm 
     1076  
     1077              IF ( sec%llstrpond ) THEN 
     1078                 transports_3d_h(2,jsec,jseg,jk) = transports_3d_h(2,jsec,jseg,jk)  + zTnorm * zrhoi 
     1079                 transports_3d_h(3,jsec,jseg,jk) = transports_3d_h(3,jsec,jseg,jk)  + zTnorm * zrhop 
     1080                 transports_3d_h(4,jsec,jseg,jk) = transports_3d_h(4,jsec,jseg,jk)  + zTnorm * 4.e+3_wp * (ztn+273.15) * 1026._wp 
     1081                 transports_3d_h(5,jsec,jseg,jk) = transports_3d_h(5,jsec,jseg,jk)  + zTnorm * 0.001 * zsn * 1026._wp 
     1082              ENDIF 
     1083 
     1084           ENDDO !end of loop on the level 
     1085 
     1086#if defined key_lim2 || defined key_lim3 
     1087 
     1088           !ICE CASE     
     1089           !------------ 
     1090           IF( sec%ll_ice_section )THEN 
     1091              SELECT CASE (sec%direction(jseg)) 
     1092              CASE(0) 
     1093                 zumid_ice = 0 
     1094                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1)) 
     1095              CASE(1) 
     1096                 zumid_ice = 0 
     1097                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1)) 
     1098              CASE(2) 
     1099                 zvmid_ice = 0 
     1100                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1)) 
     1101              CASE(3) 
     1102                 zvmid_ice = 0 
     1103                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1)) 
     1104              END SELECT 
     1105    
     1106              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J) 
     1107    
     1108              transports_2d_h(1,jsec,jseg) = transports_2d_h(1,jsec,jseg) + (zTnorm)*   & 
     1109                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
     1110                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  & 
     1111                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 
     1112                                   +zice_vol_pos 
     1113              transports_2d_h(2,jsec,jseg) = transports_2d_h(2,jsec,jseg) + (zTnorm)*   & 
     1114                                    (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
     1115                                   +zice_surf_pos 
     1116    
     1117           ENDIF !end of ice case 
     1118#endif 
     1119  
     1120        ENDDO !end of loop on the segment 
     1121 
     1122     ENDIF   !end of sec%nb_point =0 case 
     1123     ! 
     1124  END SUBROUTINE transport_h 
     1125  
    8111126  SUBROUTINE dia_dct_sum(sec,jsec)  
    8121127     !!-------------------------------------------------------------  
    … …  
    8901205              SELECT CASE( sec%direction(jseg) )  
    8911206              CASE(0,1)  
    892                  ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )  
    893                  zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )  
    894                  zrhop = interp(k%I,k%J,jk,'V',rhop)  
    895                  zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)  
     1207                 ztn   = interp(k%I,k%J,jk,'V',0 )  
     1208                 zsn   = interp(k%I,k%J,jk,'V',1 )  
     1209                 zrhop = interp(k%I,k%J,jk,'V',2)  
     1210                 zrhoi = interp(k%I,k%J,jk,'V',3)  
    8961211 
    8971212              CASE(2,3)  
    898                  ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )  
    899                  zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )  
    900                  zrhop = interp(k%I,k%J,jk,'U',rhop)  
    901                  zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)  
     1213                 ztn   = interp(k%I,k%J,jk,'U',0 )  
     1214                 zsn   = interp(k%I,k%J,jk,'U',1 )  
     1215                 zrhop = interp(k%I,k%J,jk,'U',2 )  
     1216                 zrhoi = interp(k%I,k%J,jk,'U',3 )  
    9021217                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)   
    9031218              END SELECT  
    … …  
    9571272                          sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk)  
    9581273                       ENDIF  
     1274 
     1275                       IF ( transports_3d(4,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1276                          sec%transport(7,jclass) = sec%transport(7,jclass)+transports_3d(4,jsec,jseg,jk) 
     1277                       ELSE 
     1278                          sec%transport(8,jclass) = sec%transport(8,jclass)+transports_3d(4,jsec,jseg,jk) 
     1279                       ENDIF 
     1280 
     1281                       IF ( transports_3d(5,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1282                          sec%transport( 9,jclass) = sec%transport( 9,jclass)+transports_3d(5,jsec,jseg,jk) 
     1283                       ELSE 
     1284                          sec%transport(10,jclass) = sec%transport(10,jclass)+transports_3d(5,jsec,jseg,jk) 
     1285                       ENDIF 
    9591286  
    9601287                    ELSE  
    … …  
    9631290                       sec%transport( 5,jclass) = 0._wp  
    9641291                       sec%transport( 6,jclass) = 0._wp  
     1292                       sec%transport( 7,jclass) = 0._wp 
     1293                       sec%transport( 8,jclass) = 0._wp 
     1294                       sec%transport( 9,jclass) = 0._wp 
     1295                       sec%transport(10,jclass) = 0._wp 
    9651296                    ENDIF  
    9661297  
    … …  
    9771308  
    9781309              IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN  
    979                  sec%transport( 7,1) = sec%transport( 7,1)+transports_2d(1,jsec,jseg)*1.E-6  
     1310                 sec%transport(11,1) = sec%transport(11,1)+transports_2d(1,jsec,jseg)*1.E-6  
    9801311              ELSE  
    981                  sec%transport( 8,1) = sec%transport( 8,1)+transports_2d(1,jsec,jseg)*1.E-6  
     1312                 sec%transport(12,1) = sec%transport(12,1)+transports_2d(1,jsec,jseg)*1.E-6  
    9821313              ENDIF  
    9831314  
    9841315              IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN  
    985                  sec%transport( 9,1) = sec%transport( 9,1)+transports_2d(2,jsec,jseg)*1.E-6  
     1316                 sec%transport(13,1) = sec%transport(13,1)+transports_2d(2,jsec,jseg)*1.E-6  
    9861317              ELSE  
    987                  sec%transport(10,1) = sec%transport(10,1)+transports_2d(2,jsec,jseg)*1.E-6  
     1318                 sec%transport(14,1) = sec%transport(14,1)+transports_2d(2,jsec,jseg)*1.E-6  
    9881319              ENDIF  
    9891320  
    … …  
    9951326     ELSE  !if sec%nb_point =0  
    9961327        sec%transport(1:2,:)=0.  
    997         IF (sec%llstrpond) sec%transport(3:6,:)=0.  
    998         IF (sec%ll_ice_section) sec%transport(7:10,:)=0.  
     1328        IF (sec%llstrpond) sec%transport(3:10,:)=0.  
     1329        IF (sec%ll_ice_section) sec%transport(11:14,:)=0.  
    9991330     ENDIF !end of sec%nb_point =0 case  
    10001331  
    10011332  END SUBROUTINE dia_dct_sum  
     1333 
     1334  SUBROUTINE dia_dct_sum_h(sec,jsec) 
     1335     !!------------------------------------------------------------- 
     1336     !! Exactly as dia_dct_sum but for hourly files containing data summed at each time step 
     1337     !! 
     1338     !! Purpose: Average the transport over nn_dctwri time steps  
     1339     !! and sum over the density/salinity/temperature/depth classes 
     1340     !! 
     1341     !! Method:  
     1342     !!           Sum over relevant grid cells to obtain values 
     1343     !!           for each 
     1344     !!              There are several loops:                  
     1345     !!              loop on the segment between 2 nodes 
     1346     !!              loop on the level jk 
     1347     !!              loop on the density/temperature/salinity/level classes 
     1348     !!              test on the density/temperature/salinity/level 
     1349     !! 
     1350     !!  ** Method  :Transport through a given section is equal to the sum of transports 
     1351     !!              computed on each proc. 
     1352     !!              On each proc,transport is equal to the sum of transport computed through 
     1353     !!              segments linking each point of sec%listPoint  with the next one.    
     1354     !! 
     1355     !!------------------------------------------------------------- 
     1356     !! * arguments 
     1357     TYPE(SECTION),INTENT(INOUT) :: sec 
     1358     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section 
     1359 
     1360     TYPE(POINT_SECTION) :: k 
     1361     INTEGER  :: jk,jseg,jclass                        !loop on level/segment/classes  
     1362     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 
     1363     !!------------------------------------------------------------- 
     1364 
     1365     !! Sum the relevant segments to obtain values for each class 
     1366     IF(sec%nb_point .NE. 0)THEN    
     1367 
     1368        !--------------------------------------! 
     1369        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  ! 
     1370        !--------------------------------------! 
     1371        DO jseg=1,MAX(sec%nb_point-1,0) 
     1372            
     1373           !------------------------------------------------------------------------------------------- 
     1374           ! Select the appropriate coordinate for computing the velocity of the segment 
     1375           ! 
     1376           !                      CASE(0)                                    Case (2) 
     1377           !                      -------                                    -------- 
     1378           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)       
     1379           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               | 
     1380           !                                                                            | 
     1381           !                                                                            | 
     1382           !                                                                            | 
     1383           !                      Case (3)                                            U(i,j) 
     1384           !                      --------                                              | 
     1385           !                                                                            | 
     1386           !  listPoint(jseg+1) F(i,j+1)                                                | 
     1387           !                        |                                                   | 
     1388           !                        |                                                   | 
     1389           !                        |                                 listPoint(jseg+1) F(i,j-1) 
     1390           !                        |                                             
     1391           !                        |                                             
     1392           !                     U(i,j+1)                                             
     1393           !                        |                                       Case(1)      
     1394           !                        |                                       ------       
     1395           !                        |                                             
     1396           !                        |                 listPoint(jseg+1)             listPoint(jseg)                            
     1397           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                            
     1398           ! listPoint(jseg)     F(i,j) 
     1399           !  
     1400           !------------------------------------------------------------------------------------------- 
     1401 
     1402           SELECT CASE( sec%direction(jseg) ) 
     1403           CASE(0)  ;   k = sec%listPoint(jseg) 
     1404           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 
     1405           CASE(2)  ;   k = sec%listPoint(jseg) 
     1406           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 
     1407           END SELECT 
     1408 
     1409           !---------------------------| 
     1410           !     LOOP ON THE LEVEL     | 
     1411           !---------------------------| 
     1412           !Sum of the transport on the vertical  
     1413           DO jk=1,mbathy(k%I,k%J) 
     1414 
     1415              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point 
     1416              SELECT CASE( sec%direction(jseg) ) 
     1417              CASE(0,1) 
     1418                 ztn   = interp(k%I,k%J,jk,'V',0 ) 
     1419                 zsn   = interp(k%I,k%J,jk,'V',1 ) 
     1420                 zrhop = interp(k%I,k%J,jk,'V',2 ) 
     1421                 zrhoi = interp(k%I,k%J,jk,'V',3 ) 
     1422                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
     1423              CASE(2,3) 
     1424                 ztn   = interp(k%I,k%J,jk,'U',0 ) 
     1425                 zsn   = interp(k%I,k%J,jk,'U',1 ) 
     1426                 zrhop = interp(k%I,k%J,jk,'U',2 ) 
     1427                 zrhoi = interp(k%I,k%J,jk,'U',3 ) 
     1428                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
     1429              END SELECT 
     1430 
     1431              zfsdep= fsdept(k%I,k%J,jk) 
     1432  
     1433              !------------------------------- 
     1434              !  LOOP ON THE DENSITY CLASSES | 
     1435              !------------------------------- 
     1436              !The computation is made for each density/heat/salt/... class 
     1437              DO jclass=1,MAX(1,sec%nb_class-1) 
     1438 
     1439                 !----------------------------------------------! 
     1440                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!  
     1441                 !----------------------------------------------! 
     1442  
     1443                 IF ( (                                                    & 
     1444                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      & 
     1445                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     & 
     1446                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   & 
     1447 
     1448                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    & 
     1449                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     & 
     1450                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   & 
     1451 
     1452                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   & 
     1453                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 & 
     1454                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    & 
     1455 
     1456                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   & 
     1457                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 & 
     1458                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     & 
     1459 
     1460                    ((( zfsdep .GE. sec%zlay(jclass)) .AND.                & 
     1461                    (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              & 
     1462                    ( sec%zlay(jclass) .EQ. 99. ))                         & 
     1463                                                                   ))   THEN 
     1464 
     1465                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS 
     1466                    !---------------------------------------------------------------------------- 
     1467                    IF (transports_3d_h(1,jsec,jseg,jk) .GE. 0.0) THEN  
     1468                       sec%transport_h(1,jclass) = sec%transport_h(1,jclass)+transports_3d_h(1,jsec,jseg,jk) 
     1469                    ELSE 
     1470                       sec%transport_h(2,jclass) = sec%transport_h(2,jclass)+transports_3d_h(1,jsec,jseg,jk) 
     1471                    ENDIF 
     1472                    IF( sec%llstrpond )THEN 
     1473 
     1474                       IF ( transports_3d_h(2,jsec,jseg,jk) .GE. 0.0 ) THEN  
     1475                          sec%transport_h(3,jclass) = sec%transport_h(3,jclass)+transports_3d_h(2,jsec,jseg,jk)  
     1476                       ELSE  
     1477                          sec%transport_h(4,jclass) = sec%transport_h(4,jclass)+transports_3d_h(2,jsec,jseg,jk)  
     1478                       ENDIF  
     1479  
     1480                       IF ( transports_3d_h(3,jsec,jseg,jk) .GE. 0.0 ) THEN  
     1481                          sec%transport_h(5,jclass) = sec%transport_h(5,jclass)+transports_3d_h(3,jsec,jseg,jk)  
     1482                       ELSE  
     1483                          sec%transport_h(6,jclass) = sec%transport_h(6,jclass)+transports_3d_h(3,jsec,jseg,jk)  
     1484                       ENDIF  
     1485 
     1486                       IF ( transports_3d_h(4,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1487                          sec%transport_h(7,jclass) = sec%transport_h(7,jclass)+transports_3d_h(4,jsec,jseg,jk) 
     1488                       ELSE 
     1489                          sec%transport_h(8,jclass) = sec%transport_h(8,jclass)+transports_3d_h(4,jsec,jseg,jk) 
     1490                       ENDIF 
     1491 
     1492                       IF ( transports_3d_h(5,jsec,jseg,jk) .GE. 0.0 ) THEN 
     1493                          sec%transport_h( 9,jclass) = sec%transport_h( 9,jclass)+transports_3d_h(5,jsec,jseg,jk) 
     1494                       ELSE 
     1495                          sec%transport_h(10,jclass) = sec%transport_h(10,jclass)+transports_3d_h(5,jsec,jseg,jk) 
     1496                       ENDIF 
     1497 
     1498                    ELSE 
     1499                       sec%transport_h( 3,jclass) = 0._wp 
     1500                       sec%transport_h( 4,jclass) = 0._wp 
     1501                       sec%transport_h( 5,jclass) = 0._wp 
     1502                       sec%transport_h( 6,jclass) = 0._wp 
     1503                       sec%transport_h( 7,jclass) = 0._wp 
     1504                       sec%transport_h( 8,jclass) = 0._wp 
     1505                       sec%transport_h( 9,jclass) = 0._wp 
     1506                       sec%transport_h(10,jclass) = 0._wp 
     1507                    ENDIF 
     1508 
     1509                 ENDIF ! end of test if point is in class 
     1510    
     1511              ENDDO ! end of loop on the classes 
     1512 
     1513           ENDDO ! loop over jk 
     1514 
     1515#if defined key_lim2 || defined key_lim3 
     1516 
     1517           !ICE CASE     
     1518           IF( sec%ll_ice_section )THEN 
     1519 
     1520              IF ( transports_2d_h(1,jsec,jseg) .GE. 0.0 ) THEN 
     1521                 sec%transport_h(11,1) = sec%transport_h(11,1)+transports_2d_h(1,jsec,jseg) 
     1522              ELSE 
     1523                 sec%transport_h(12,1) = sec%transport_h(12,1)+transports_2d_h(1,jsec,jseg) 
     1524              ENDIF 
     1525 
     1526              IF ( transports_2d_h(3,jsec,jseg) .GE. 0.0 ) THEN 
     1527                 sec%transport_h(13,1) = sec%transport_h(13,1)+transports_2d_h(2,jsec,jseg) 
     1528              ELSE 
     1529                 sec%transport_h(14,1) = sec%transport_h(14,1)+transports_2d_h(2,jsec,jseg) 
     1530              ENDIF 
     1531 
     1532           ENDIF !end of ice case 
     1533#endif 
     1534  
     1535        ENDDO !end of loop on the segment 
     1536 
     1537     ELSE  !if sec%nb_point =0 
     1538        sec%transport_h(1:2,:)=0. 
     1539        IF (sec%llstrpond) sec%transport_h(3:10,:)=0. 
     1540        IF (sec%ll_ice_section) sec%transport_h( 11:14,:)=0. 
     1541     ENDIF !end of sec%nb_point =0 case 
     1542 
     1543  END SUBROUTINE dia_dct_sum_h 
    10021544   
     1545  SUBROUTINE dia_dct_wri_NOOS(kt,ksec,sec) 
     1546     !!------------------------------------------------------------- 
     1547     !! Write transport output in numdct using NOOS formatting  
     1548     !!  
     1549     !! Purpose: Write  transports in ascii files 
     1550     !!  
     1551     !! Method: 
     1552     !!        1. Write volume transports in "volume_transport" 
     1553     !!           Unit: Sv : area * Velocity / 1.e6  
     1554     !!  
     1555     !!        2. Write heat transports in "heat_transport" 
     1556     !!           Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15 
     1557     !!  
     1558     !!        3. Write salt transports in "salt_transport" 
     1559     !!           Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6 
     1560     !! 
     1561     !!-------------------------------------------------------------  
     1562     !!arguments 
     1563     INTEGER, INTENT(IN)          :: kt          ! time-step 
     1564     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write    
     1565     INTEGER ,INTENT(IN)          :: ksec        ! section number 
     1566 
     1567     !!local declarations 
     1568     INTEGER               :: jclass,ji             ! Dummy loop 
     1569     CHARACTER(len=2)      :: classe             ! Classname  
     1570     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds 
     1571     REAL(wp)              :: zslope             ! section's slope coeff 
     1572     ! 
     1573     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace  
     1574     !!-------------------------------------------------------------  
     1575     CALL wrk_alloc(nb_type_class , zsumclasses )   
     1576 
     1577     zsumclasses(:)=0._wp 
     1578     zslope = sec%slopeSection        
     1579 
     1580     WRITE(numdct_NOOS,'(I4,a1,I2,a1,I2,a12,i3,a17,i3,a10,a25)') nyear,'.',nmonth,'.',nday,'   Transect:',ksec-1,'   No. of layers:',sec%nb_class-1,'   Name: ',sec%name 
     1581 
     1582     DO jclass=1,MAX(1,sec%nb_class-1) 
     1583        zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass) 
     1584     ENDDO 
     1585  
     1586     classe   = 'total   ' 
     1587     zbnd1   = 0._wp 
     1588     zbnd2   = 0._wp 
     1589 
     1590     IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 
     1591        WRITE(numdct_NOOS,'(9e12.4E2)') -(zsumclasses( 1)+zsumclasses( 2)), -zsumclasses( 2),-zsumclasses( 1),   & 
     1592                                        -(zsumclasses( 7)+zsumclasses( 8)), -zsumclasses( 8),-zsumclasses( 7),   & 
     1593                                        -(zsumclasses( 9)+zsumclasses(10)), -zsumclasses(10),-zsumclasses( 9) 
     1594     ELSE 
     1595        WRITE(numdct_NOOS,'(9e12.4E2)')   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   & 
     1596                                          zsumclasses( 7)+zsumclasses( 8) ,  zsumclasses( 7), zsumclasses( 8),   & 
     1597                                          zsumclasses( 9)+zsumclasses(10) ,  zsumclasses( 9), zsumclasses(10) 
     1598     ENDIF  
     1599 
     1600     DO jclass=1,MAX(1,sec%nb_class-1) 
     1601    
     1602        classe   = 'N       ' 
     1603        zbnd1   = 0._wp 
     1604        zbnd2   = 0._wp 
     1605 
     1606        !insitu density classes transports 
     1607        IF( ( sec%zsigi(jclass)   .NE. 99._wp ) .AND. & 
     1608            ( sec%zsigi(jclass+1) .NE. 99._wp )       )THEN 
     1609           classe = 'DI       ' 
     1610           zbnd1 = sec%zsigi(jclass) 
     1611           zbnd2 = sec%zsigi(jclass+1) 
     1612        ENDIF 
     1613        !potential density classes transports 
     1614        IF( ( sec%zsigp(jclass)   .NE. 99._wp ) .AND. & 
     1615            ( sec%zsigp(jclass+1) .NE. 99._wp )       )THEN 
     1616           classe = 'DP      ' 
     1617           zbnd1 = sec%zsigp(jclass) 
     1618           zbnd2 = sec%zsigp(jclass+1) 
     1619        ENDIF 
     1620        !depth classes transports 
     1621        IF( ( sec%zlay(jclass)    .NE. 99._wp ) .AND. & 
     1622            ( sec%zlay(jclass+1)  .NE. 99._wp )       )THEN  
     1623           classe = 'Z       ' 
     1624           zbnd1 = sec%zlay(jclass) 
     1625           zbnd2 = sec%zlay(jclass+1) 
     1626        ENDIF 
     1627        !salinity classes transports 
     1628        IF( ( sec%zsal(jclass) .NE. 99._wp    ) .AND. & 
     1629            ( sec%zsal(jclass+1) .NE. 99._wp  )       )THEN 
     1630           classe = 'S       ' 
     1631           zbnd1 = sec%zsal(jclass) 
     1632           zbnd2 = sec%zsal(jclass+1)    
     1633        ENDIF 
     1634        !temperature classes transports 
     1635        IF( ( sec%ztem(jclass) .NE. 99._wp     ) .AND. & 
     1636            ( sec%ztem(jclass+1) .NE. 99._wp     )       ) THEN 
     1637           classe = 'T       ' 
     1638           zbnd1 = sec%ztem(jclass) 
     1639           zbnd2 = sec%ztem(jclass+1) 
     1640        ENDIF 
     1641                   
     1642        !write volume transport per class 
     1643        IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 
     1644           WRITE(numdct_NOOS,'(9e12.4E2)') -(sec%transport( 1,jclass)+sec%transport( 2,jclass)),-sec%transport( 2,jclass),-sec%transport( 1,jclass), & 
     1645                                           -(sec%transport( 7,jclass)+sec%transport( 8,jclass)),-sec%transport( 8,jclass),-sec%transport( 7,jclass), & 
     1646                                           -(sec%transport( 9,jclass)+sec%transport(10,jclass)),-sec%transport(10,jclass),-sec%transport( 9,jclass) 
     1647        ELSE 
     1648           WRITE(numdct_NOOS,'(9e12.4E2)')   sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), & 
     1649                                             sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), & 
     1650                                             sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass) 
     1651        ENDIF 
     1652 
     1653     ENDDO 
     1654 
     1655     CALL wrk_dealloc(nb_type_class , zsumclasses )   
     1656 
     1657  END SUBROUTINE dia_dct_wri_NOOS 
     1658 
     1659  SUBROUTINE dia_dct_wri_NOOS_h(hr,ksec,sec) 
     1660     !!------------------------------------------------------------- 
     1661     !! As routine dia_dct_wri_NOOS but for hourly output files 
     1662     !! 
     1663     !! Write transport output in numdct using NOOS formatting  
     1664     !!  
     1665     !! Purpose: Write  transports in ascii files 
     1666     !!  
     1667     !! Method: 
     1668     !!        1. Write volume transports in "volume_transport" 
     1669     !!           Unit: Sv : area * Velocity / 1.e6  
     1670     !! 
     1671     !!-------------------------------------------------------------  
     1672     !!arguments 
     1673     INTEGER, INTENT(IN)          :: hr          ! hour 
     1674     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write    
     1675     INTEGER ,INTENT(IN)          :: ksec        ! section number 
     1676 
     1677     !!local declarations 
     1678     INTEGER               :: jclass,jhr            ! Dummy loop 
     1679     CHARACTER(len=2)      :: classe             ! Classname  
     1680     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds 
     1681     REAL(wp)              :: zslope             ! section's slope coeff 
     1682     ! 
     1683     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace  
     1684     !!-------------------------------------------------------------  
     1685 
     1686     CALL wrk_alloc(nb_type_class , zsumclasses )  
     1687 
     1688     zsumclasses(:)=0._wp 
     1689     zslope = sec%slopeSection        
     1690 
     1691     DO jclass=1,MAX(1,sec%nb_class-1) 
     1692        zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport_h(1:nb_type_class,jclass) 
     1693     ENDDO 
     1694  
     1695     !write volume transport per class 
     1696     IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 
     1697        z_hr_output(ksec,hr,1)=-(zsumclasses(1)+zsumclasses(2)) 
     1698     ELSE 
     1699        z_hr_output(ksec,hr,1)= (zsumclasses(1)+zsumclasses(2)) 
     1700     ENDIF 
     1701 
     1702     DO jclass=1,MAX(1,sec%nb_class-1) 
     1703        IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 
     1704           z_hr_output(ksec,hr,jclass+1)=-(sec%transport_h(1,jclass)+sec%transport_h(2,jclass)) 
     1705        ELSE 
     1706           z_hr_output(ksec,hr,jclass+1)= (sec%transport_h(1,jclass)+sec%transport_h(2,jclass)) 
     1707        ENDIF 
     1708     ENDDO 
     1709 
     1710     IF ( hr .eq. 48._wp ) THEN 
     1711        WRITE(numdct_NOOS_h,'(I4,a1,I2,a1,I2,a12,i3,a17,i3)') nyear,'.',nmonth,'.',nday,'   Transect:',ksec-1,'   No. of layers:',sec%nb_class-1 
     1712        DO jhr=25,48 
     1713           WRITE(numdct_NOOS_h,'(11F12.1)')  z_hr_output(ksec,jhr,1), (z_hr_output(ksec,jhr,jclass+1), jclass=1,MAX(1,10) ) 
     1714        ENDDO 
     1715     ENDIF 
     1716 
     1717     CALL wrk_dealloc(nb_type_class , zsumclasses ) 
     1718 
     1719  END SUBROUTINE dia_dct_wri_NOOS_h 
     1720 
    10031721  SUBROUTINE dia_dct_wri(kt,ksec,sec) 
    10041722     !!------------------------------------------------------------- 
    … …  
    10921810           WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope,  & 
    10931811                              jclass,classe,zbnd1,zbnd2,& 
    1094                               sec%transport(3,jclass)*1.e-15,sec%transport(4,jclass)*1.e-15, & 
    1095                               ( sec%transport(3,jclass)+sec%transport(4,jclass) )*1.e-15 
     1812                              sec%transport(7,jclass)*1.e-15,sec%transport(8,jclass)*1.e-15, & 
     1813                              ( sec%transport(7,jclass)+sec%transport(8,jclass) )*1.e-15 
    10961814           !write salt transport per class 
    10971815           WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope,  & 
    10981816                              jclass,classe,zbnd1,zbnd2,& 
    1099                               sec%transport(5,jclass)*1.e-9,sec%transport(6,jclass)*1.e-9,& 
    1100                               (sec%transport(5,jclass)+sec%transport(6,jclass))*1.e-9 
     1817                              sec%transport(9,jclass)*1.e-9,sec%transport(10,jclass)*1.e-9,& 
     1818                              (sec%transport(9,jclass)+sec%transport(10,jclass))*1.e-9 
    11011819        ENDIF 
    11021820 
    … …  
    11171835        WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, & 
    11181836                           jclass,"total",zbnd1,zbnd2,& 
    1119                            zsumclasses(3)*1.e-15,zsumclasses(4)*1.e-15,& 
    1120                            (zsumclasses(3)+zsumclasses(4) )*1.e-15 
     1837                           zsumclasses(7)*1.e-15,zsumclasses(8)*1.e-15,& 
     1838                           (zsumclasses(7)+zsumclasses(8) )*1.e-15 
    11211839        !write total salt transport 
    11221840        WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, & 
    11231841                           jclass,"total",zbnd1,zbnd2,& 
    1124                            zsumclasses(5)*1.e-9,zsumclasses(6)*1.e-9,& 
    1125                            (zsumclasses(5)+zsumclasses(6))*1.e-9 
     1842                           zsumclasses(9)*1.e-9,zsumclasses(10)*1.e-9,& 
     1843                           (zsumclasses(9)+zsumclasses(10))*1.e-9 
    11261844     ENDIF 
    11271845 
    … …  
    11311849        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,& 
    11321850                              jclass,"ice_vol",zbnd1,zbnd2,& 
    1133                               sec%transport(7,1),sec%transport(8,1),& 
    1134                               sec%transport(7,1)+sec%transport(8,1) 
     1851                              sec%transport(11,1),sec%transport(12,1),& 
     1852                              sec%transport(11,1)+sec%transport(12,1) 
    11351853        !write total ice surface transport 
    11361854        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,& 
    11371855                              jclass,"ice_surf",zbnd1,zbnd2,& 
    1138                               sec%transport(9,1),sec%transport(10,1), & 
    1139                               sec%transport(9,1)+sec%transport(10,1)  
     1856                              sec%transport(13,1),sec%transport(14,1), & 
     1857                              sec%transport(13,1)+sec%transport(14,1)  
    11401858     ENDIF 
    11411859                                               
    … …  
    11461864  END SUBROUTINE dia_dct_wri 
    11471865 
    1148   FUNCTION interp(ki, kj, kk, cd_point, ptab) 
     1866   PURE  FUNCTION interp(ki, kj, kk, cd_point,var)  
    11491867  !!---------------------------------------------------------------------- 
    11501868  !! 
    … …  
    12081926  !*arguments 
    12091927  INTEGER, INTENT(IN)                          :: ki, kj, kk   ! coordinate of point 
     1928  INTEGER, INTENT(IN)                          :: var   !  which variable 
    12101929  CHARACTER(len=1), INTENT(IN)                 :: cd_point     ! type of point (U, V) 
    1211   REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab         ! variable to compute at (ki, kj, kk ) 
    12121930  REAL(wp)                                     :: interp       ! interpolated variable  
    12131931 
    … …  
    12201938  !!---------------------------------------------------------------------- 
    12211939 
     1940  
     1941 
    12221942  IF( cd_point=='U' )THEN  
    12231943     ii1 = ki    ; ij1 = kj  
    … …  
    12501970   
    12511971     ! result 
    1252      interp = zmsk * ( zwgt2 *  ptab(ii1,ij1,kk) + zwgt1 *  ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )    
    1253  
     1972           SELECT CASE( var ) 
     1973              CASE(0)  ;    interp = zmsk * ( zwgt2 *  tsn(ii1,ij1,kk,jp_tem) + zwgt1 *  tsn(ii1,ij1,kk,jp_tem) ) / ( zwgt2 + zwgt1 ) 
     1974              CASE(1)  ;    interp = zmsk * ( zwgt2 *  tsn(ii1,ij1,kk,jp_sal) + zwgt1 *  tsn(ii1,ij1,kk,jp_sal) ) / ( zwgt2 + zwgt1 ) 
     1975              CASE(2)  ;    interp = zmsk * ( zwgt2 *  rhop(ii1,ij1,kk) + zwgt1 *  rhop(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 ) 
     1976              CASE(3)  ;    interp = zmsk * ( zwgt2 *  (rhd(ii1,ij1,kk)*rau0+rau0) + zwgt1 *  (rhd(ii1,ij1,kk)*rau0+rau0) ) / ( zwgt2 + zwgt1 ) 
     1977           END SELECT 
    12541978 
    12551979  ELSE       ! full step or partial step case  
    … …  
    12731997        IF( ze3t >= 0. )THEN  
    12741998           ! zbis 
    1275            zbis = ptab(ii2,ij2,kk) + zwgt1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) )  
     1999           SELECT CASE( var ) 
     2000           CASE(0)   
     2001                     zbis   =  tsn(ii2,ij2,kk,jp_tem) + zwgt1 *  (tsn(ii2,ij2,kk-1,jp_tem)-tsn(ii2,ij2,kk,jp_tem)   ) 
     2002                     interp =  zmsk * ( zet2 * tsn(ii1,ij1,kk,jp_tem) + zet1 * zbis )/( zet1 + zet2 ) 
     2003           CASE(1)   
     2004                     zbis   =  tsn(ii2,ij2,kk,jp_sal) + zwgt1 *  (tsn(ii2,ij2,kk-1,jp_sal)-tsn(ii2,ij2,kk,jp_sal)   ) 
     2005                     interp =  zmsk * ( zet2 * tsn(ii1,ij1,kk,jp_sal) + zet1 * zbis )/( zet1 + zet2 ) 
     2006           CASE(2)   
     2007                     zbis   =  rhop(ii2,ij2,kk) + zwgt1 *  (rhop(ii2,ij2,kk-1)-rhop(ii2,ij2,kk)   ) 
     2008                     interp =  zmsk * ( zet2 * rhop(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 ) 
     2009           CASE(3)   
     2010                     zbis   =  (rhd(ii2,ij2,kk)*rau0+rau0) + zwgt1 *  ( (rhd(ii2,ij2,kk-1)*rau0+rau0)-(rhd(ii2,ij2,kk)*rau0+rau0)   ) 
     2011                     interp =  zmsk * ( zet2 * (rhd(ii1,ij1,kk)*rau0+rau0) + zet1 * zbis )/( zet1 + zet2 ) 
     2012           END SELECT 
    12762013           ! result 
    1277             interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 ) 
    12782014        ELSE 
    12792015           ! zbis 
    1280            zbis = ptab(ii1,ij1,kk) + zwgt2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) ) 
    1281            ! result 
    1282            interp = zmsk * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 ) 
     2016           SELECT CASE( var ) 
     2017           CASE(0)   
     2018                 zbis   = tsn(ii1,ij1,kk,jp_tem) + zwgt2 * ( tsn(ii1,ij1,kk-1,jp_tem) - tsn(ii1,ij2,kk,jp_tem) ) 
     2019                 interp = zmsk * ( zet2 * zbis + zet1 * tsn(ii2,ij2,kk,jp_tem) )/( zet1 + zet2 ) 
     2020           CASE(1)   
     2021                 zbis   = tsn(ii1,ij1,kk,jp_sal) + zwgt2 * ( tsn(ii1,ij1,kk-1,jp_sal) - tsn(ii1,ij2,kk,jp_sal) ) 
     2022                 interp = zmsk * ( zet2 * zbis + zet1 * tsn(ii2,ij2,kk,jp_sal) )/( zet1 + zet2 ) 
     2023           CASE(2)   
     2024                 zbis   = rhop(ii1,ij1,kk) + zwgt2 * ( rhop(ii1,ij1,kk-1) - rhop(ii1,ij2,kk) ) 
     2025                 interp = zmsk * ( zet2 * zbis + zet1 * rhop(ii2,ij2,kk) )/( zet1 + zet2 ) 
     2026           CASE(3)   
     2027                 zbis   = (rhd(ii1,ij1,kk)*rau0+rau0) + zwgt2 * ( (rhd(ii1,ij1,kk-1)*rau0+rau0) - (rhd(ii1,ij2,kk)*rau0+rau0) ) 
     2028                 interp = zmsk * ( zet2 * zbis + zet1 * (rhd(ii2,ij2,kk)*rau0+rau0) )/( zet1 + zet2 ) 
     2029           END SELECT 
    12832030        ENDIF     
    12842031 
    12852032     ELSE 
    1286         interp = zmsk * (  zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 ) 
     2033        SELECT CASE( var ) 
     2034        CASE(0)   
     2035           interp = zmsk * (  zet2 * tsn(ii1,ij1,kk,jp_tem) + zet1 * tsn(ii2,ij2,kk,jp_tem) )/( zet1 + zet2 ) 
     2036        CASE(1)   
     2037           interp = zmsk * (  zet2 * tsn(ii1,ij1,kk,jp_sal) + zet1 * tsn(ii2,ij2,kk,jp_sal) )/( zet1 + zet2 ) 
     2038        CASE(2)   
     2039           interp = zmsk * (  zet2 * rhop(ii1,ij1,kk) + zet1 * rhop(ii2,ij2,kk) )/( zet1 + zet2 ) 
     2040        CASE(3)   
     2041           interp = zmsk * (  zet2 * (rhd(ii1,ij1,kk)*rau0+rau0) + zet1 * (rhd(ii2,ij2,kk)*rau0+rau0) )/( zet1 + zet2 ) 
     2042        END SELECT 
    12872043     ENDIF 
    12882044 
    12892045  ENDIF 
    1290  
    12912046 
    12922047  END FUNCTION interp 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r10248 r10249  
    4444   USE in_out_manager  ! I/O manager 
    4545   USE diadimg         ! dimg direct access file format output 
     46   USE diatmb          ! Top,middle,bottom output 
     47   USE dia25h          ! 25h Mean output 
     48   USE diaopfoam       ! Diaopfoam output 
    4649   USE iom 
    4750   USE ioipsl 
    4851   USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities      
     52   USE eosbn2         ! equation of state                (eos_bn2 routine) 
     53 
    4954 
    5055#if defined key_lim2 
    … …  
    132137      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace 
    133138      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace 
     139      REAL(wp), POINTER, DIMENSION(:,:,:) :: zrhd , zrhop  ! 3D workspace 
    134140      !!---------------------------------------------------------------------- 
    135141      !  
    … …  
    138144      CALL wrk_alloc( jpi , jpj      , z2d ) 
    139145      CALL wrk_alloc( jpi , jpj, jpk , z3d ) 
     146      CALL wrk_alloc( jpi , jpj, jpk , zrhd , zrhop ) 
    140147      ! 
    141148      ! Output the initial state and forcings 
    … …  
    376383         CALL iom_put( "v_salttr", 0.5 * z2d )            !  heat transport in j-direction 
    377384      ENDIF 
     385 
     386      IF( iom_use("rhop") ) THEN 
     387         CALL eos( tsn, zrhd, zrhop, fsdept_n(:,:,:) )       ! now in situ and potential density 
     388         zrhop(:,:,jpk) = 0._wp 
     389         CALL iom_put( 'rhop', zrhop ) 
     390      ENDIF 
     391 
    378392      ! 
    379393      CALL wrk_dealloc( jpi , jpj      , z2d ) 
    380394      CALL wrk_dealloc( jpi , jpj, jpk , z3d ) 
     395      CALL wrk_dealloc( jpi , jpj, jpk , zrhd , zrhop ) 
     396      ! 
     397      ! If we want tmb values  
     398 
     399      IF (ln_diatmb) THEN 
     400         CALL dia_tmb 
     401      ENDIF 
     402      IF (ln_dia25h) THEN 
     403         CALL dia_25h( kt ) 
     404      ENDIF 
     405      IF (ln_diaopfoam) THEN 
     406         CALL dia_diaopfoam 
     407      ENDIF 
    381408      ! 
    382409      IF( nn_timing == 1 )   CALL timing_stop('dia_wri') 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r10246 r10249  
    135135      !!---------------------------------------------------------------------- 
    136136      USE ioipsl 
     137      NAMELIST/namrun/ ln_NOOS 
    137138      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               & 
    138          &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   & 
     139         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstdate, ln_rstart , nn_rstctl,   & 
    139140         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   & 
    140141         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler 
    … …  
    174175         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir 
    175176         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart 
     177         WRITE(numout,*) '      datestamping of restarts        ln_rstdate  = ', ln_rstdate 
    176178         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler 
    177179         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl 
    … …  
    192194         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber 
    193195         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz 
     196         WRITE(numout,*) '      NOOS transect diagnostics       ln_NOOS    = ', ln_NOOS 
    194197      ENDIF 
    195198 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r10248 r10249  
    3131   USE wrk_nemo        ! Memory allocation 
    3232   USE timing          ! Timing 
     33   USE iom    ! slwa 
    3334 
    3435   IMPLICIT NONE 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90

    r10246 r10249  
    209209      ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:)  
    210210      ! 
    211       IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==! 
     211      IF( ln_sco .and. 1==0) THEN                   !==   s- or mixed s-zps-coordinate   ==! 
    212212         ! 
    213213         CALL wrk_alloc( jpk, ztp, zsp ) 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90

    r10246 r10249  
    2424   USE wrk_nemo        ! Memory Allocation 
    2525   USE timing          ! Timing 
     26#if defined key_bdy  
     27   USE bdy_oce         ! ocean open boundary conditions  
     28#endif  
    2629 
    2730   IMPLICIT NONE 
    … …  
    7881      REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke 
    7982      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv  
     83#if defined key_bdy  
     84      INTEGER  ::   jb                 ! dummy loop indices  
     85      INTEGER  ::   ii, ij, igrd, ib_bdy   ! local integers  
     86      INTEGER  ::   fu, fv  
     87#endif  
    8088      !!---------------------------------------------------------------------- 
    8189      ! 
    … …  
    97105       
    98106      zhke(:,:,jpk) = 0._wp 
    99        
     107 
     108#if defined key_bdy  
     109      ! Maria Luneva & Fred Wobus: July-2016  
     110      ! compensate for lack of turbulent kinetic energy on liquid bdy points  
     111      DO ib_bdy = 1, nb_bdy  
     112         IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN  
     113            igrd = 2           ! Copying normal velocity into points outside bdy  
     114            DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)  
     115               DO jk = 1, jpkm1  
     116                  ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)  
     117                  ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)  
     118                  fu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) )  
     119                  un(ii-fu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk)  
     120               END DO  
     121            END DO  
     122            !  
     123            igrd = 3           ! Copying normal velocity into points outside bdy  
     124            DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)  
     125               DO jk = 1, jpkm1  
     126                  ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)  
     127                  ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)  
     128                  fv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) )  
     129                  vn(ii,ij-fv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk)  
     130               END DO  
     131            END DO  
     132         ENDIF  
     133      ENDDO   
     134#endif        
     135 
    100136      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==! 
    101137      ! 
    … …  
    112148            END DO 
    113149         END DO 
     150         ! 
     151         CALL lbc_lnk( zhke, 'T', 1. ) 
    114152         ! 
    115153      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --! 
    … …  
    134172      END SELECT 
    135173      ! 
     174 
     175#if defined key_bdy  
     176      ! restore velocity masks at points outside boundary  
     177      un(:,:,:) = un(:,:,:) * umask(:,:,:)  
     178      vn(:,:,:) = vn(:,:,:) * vmask(:,:,:)  
     179#endif   
     180 
    136181      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==! 
    137182         DO jj = 2, jpjm1 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r10248 r10249  
    4141   USE timing          ! Timing     
    4242   USE sbcapr          ! surface boundary condition: atmospheric pressure 
     43   USE diatmb          ! Top,middle,bottom output 
    4344   USE dynadv, ONLY: ln_dynadv_vec 
    4445#if defined key_agrif 
    4546   USE agrif_opa_interp ! agrif 
    46 #endif 
    47 #if defined key_asminc    
    48    USE asminc          ! Assimilation increment 
    4947#endif 
    5048 
    … …  
    144142      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices 
    145143      INTEGER  ::   ikbu, ikbv, noffset      ! local integers 
     144      REAL(wp) ::   zmdi 
    146145      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars 
    147146      REAL(wp) ::   zx1, zy1, zx2, zy2         !   -      - 
    … …  
    169168      CALL wrk_alloc( jpi, jpj, zhf ) 
    170169      ! 
     170      zmdi=1.e+20                               !  missing data indicator for masking 
    171171      !                                         !* Local constant initialization 
    172172      z1_12 = 1._wp / 12._wp  
    … …  
    459459                &                        + fwfisf(:,:) + fwfisf_b(:,:)                     ) 
    460460      ENDIF 
    461 #if defined key_asminc 
    462       !                                         ! Include the IAU weighted SSH increment 
    463       IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 
    464          zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:) 
    465       ENDIF 
    466 #endif 
    467461      !                                   !* Fill boundary data arrays for AGRIF 
    468462      !                                   ! ------------------------------------ 
    … …  
    926920      CALL wrk_dealloc( jpi, jpj, zhf ) 
    927921      ! 
     922      IF ( ln_diatmb ) THEN 
     923         CALL iom_put( "baro_u" , un_b*umask(:,:,1)+zmdi*(1-umask(:,:,1 ) ) )  ! Barotropic  U Velocity 
     924         CALL iom_put( "baro_v" , vn_b*vmask(:,:,1)+zmdi*(1-vmask(:,:,1 ) ) )  ! Barotropic  V Velocity 
     925      ENDIF 
    928926      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts') 
    929927      ! 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90

    r10248 r10249  
    7474      INTEGER, INTENT(in) ::   kt                      ! time step 
    7575      !  
    76       INTEGER             ::   jk                      ! dummy loop indice 
     76      INTEGER             ::   ji, jj, jk              ! dummy loop indices 
    7777      REAL(wp)            ::   z2dt, z1_rau0           ! local scalars 
    7878      !!---------------------------------------------------------------------- 
    … …  
    9494      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog) 
    9595      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt 
     96 
     97#if defined key_asminc 
     98      !                                                ! Include the IAU weighted SSH increment 
     99      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 
     100         CALL ssh_asm_inc( kt ) 
     101#if defined key_vvl 
     102! Don't directly adjust ssh but change hdivn at all levels instead 
     103! In trasbc also add in the heat and salt content associated with these changes at each level   
     104        DO jk = 1, jpkm1                                  
     105                 hdivn(:,:,jk) = hdivn(:,:,jk) - ( ssh_iau(:,:) / ( ht_0(:,:) + 1.0 - ssmask(:,:) ) ) * ( e3t_0(:,:,jk) / fse3t_n(:,:,jk) ) * tmask(:,:,jk)  
     106        END DO 
     107        CALL lbc_lnk( hdivn, 'T', 1. ) ! Not sure that's necessary 
     108      ENDIF 
     109#endif 
     110#endif 
    96111 
    97112      !                                           !------------------------------! 
    … …  
    123138#endif 
    124139 
    125 #if defined key_asminc 
    126       !                                                ! Include the IAU weighted SSH increment 
    127       IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 
    128          CALL ssh_asm_inc( kt ) 
    129          ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:) 
    130       ENDIF 
    131 #endif 
    132  
    133140      !                                           !------------------------------! 
    134141      !                                           !           outputs            ! 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90

    r10246 r10249  
    1818   !!---------------------------------------------------------------------- 
    1919   USE par_oce        ! NEMO parameters 
     20   USE phycst         ! for rday 
    2021   USE dom_oce        ! NEMO domain 
    2122   USE in_out_manager ! NEMO IO routines 
     23   USE ioipsl, ONLY : ju2ymds    ! for calendar 
    2224   USE lib_mpp        ! NEMO MPI library, lk_mpp in particular 
    2325   USE netcdf         ! netcdf routines for IO 
    … …  
    231233      INTEGER ::   jn   ! dummy loop index 
    232234      INTEGER ::   ix_dim, iy_dim, ik_dim, in_dim 
    233       CHARACTER(len=256)     :: cl_path 
    234       CHARACTER(len=256)     :: cl_filename 
     235      INTEGER             ::   iyear, imonth, iday 
     236      REAL (wp)           ::   zsec 
     237      REAL (wp)           ::   zfjulday 
     238      CHARACTER(len=256)  :: cl_path 
     239      CHARACTER(len=256)  :: cl_filename 
     240      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character 
    235241      TYPE(iceberg), POINTER :: this 
    236242      TYPE(point)  , POINTER :: pt 
    … …  
    240246      cl_path = TRIM(cn_ocerst_outdir) 
    241247      IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/' 
     248      IF ( ln_rstdate ) THEN 
     249         zfjulday = fjulday + rdttra(1) / rday 
     250         IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday )   zfjulday = REAL(NINT(zfjulday),wp)   ! avoid truncation error 
     251         CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )            
     252         WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday 
     253      ELSE 
     254         IF( kt > 999999999 ) THEN   ;   WRITE(clkt, *       ) kt 
     255         ELSE                        ;   WRITE(clkt, '(i8.8)') kt 
     256         ENDIF 
     257      ENDIF 
    242258      IF( lk_mpp ) THEN 
    243          WRITE(cl_filename,'(A,"_icebergs_",I8.8,"_restart_",I4.4,".nc")') TRIM(cexper), kt, narea-1 
     259         WRITE(cl_filename,'(A,"_icebergs_",A,"_restart_",I4.4,".nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)), narea-1 
    244260      ELSE 
    245          WRITE(cl_filename,'(A,"_icebergs_",I8.8,"_restart.nc")') TRIM(cexper), kt 
     261         WRITE(cl_filename,'(A,"_icebergs_",A,"_restart.nc")') TRIM(cexper), TRIM(ADJUSTL(clkt)) 
    246262      ENDIF 
    247263      IF (nn_verbose_level >= 0) WRITE(numout,'(2a)') 'icebergs, write_restart: creating ',TRIM(cl_path)//TRIM(cl_filename) 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90

    r10246 r10249  
    3030   CHARACTER(lc) ::   cn_ocerst_outdir !: restart output directory 
    3131   LOGICAL       ::   ln_rstart        !: start from (F) rest or (T) a restart file 
     32   LOGICAL       ::   ln_rstdate       !: datestamping of restarts  
    3233   LOGICAL       ::   ln_rst_list      !: output restarts at list of times (T) or by frequency (F) 
    3334   INTEGER       ::   nn_no            !: job number 
    … …  
    4849   LOGICAL       ::   ln_clobber       !: clobber (overwrite) an existing file 
    4950   INTEGER       ::   nn_chunksz       !: chunksize (bytes) for NetCDF file (works only with iom_nf90 routines) 
     51   LOGICAL       ::   ln_NOOS = .FALSE.                !: NOOS transects  diagnostics 
    5052#if defined key_netcdf4 
    5153   !!---------------------------------------------------------------------- 
    … …  
    133135   INTEGER ::   numdct_heat     =   -1      !: logical unit for heat    transports output 
    134136   INTEGER ::   numdct_salt     =   -1      !: logical unit for salt    transports output 
     137   INTEGER ::   numdct_NOOS     =   -1      !: logical unit for NOOS    transports output 
     138   INTEGER ::   numdct_NOOS_h   =   -1      !: logical unit for NOOS hourly transports output 
    135139   INTEGER ::   numfl           =   -1      !: logical unit for floats ascii output 
    136140   INTEGER ::   numflo          =   -1      !: logical unit for floats ascii output 
     141    
     142   INTEGER ::   numdct_reg_bin  =   -1      !: logical unit for NOOS    transports output 
     143   INTEGER ::   numdct_reg_txt  =   -1      !: logical unit for NOOS hourly transports output 
    137144 
    138145   !!---------------------------------------------------------------------- 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/IOM/iom_def.F90

    r10246 r10249  
    4343   INTEGER, PARAMETER, PUBLIC ::   jp_i1    = 204      !: write INTEGER(1) 
    4444 
    45    INTEGER, PARAMETER, PUBLIC ::   jpmax_files  = 100   !: maximum number of simultaneously opened file 
     45   INTEGER, PARAMETER, PUBLIC ::   jpmax_files  = 200   !: maximum number of simultaneously opened file 
    4646   INTEGER, PARAMETER, PUBLIC ::   jpmax_vars   = 600  !: maximum number of variables in one file 
    4747   INTEGER, PARAMETER, PUBLIC ::   jpmax_dims   =  4   !: maximum number of dimensions for one variable 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r10246 r10249  
    2121   USE in_out_manager  ! I/O manager 
    2222   USE iom             ! I/O module 
     23   USE ioipsl, ONLY : ju2ymds    ! for calendar 
    2324   USE eosbn2          ! equation of state            (eos bn2 routine) 
    2425   USE trdmxl_oce      ! ocean active mixed layer tracers trends variables 
    … …  
    5455      !!---------------------------------------------------------------------- 
    5556      INTEGER, INTENT(in) ::   kt     ! ocean time-step 
     57      INTEGER             ::   iyear, imonth, iday 
     58      REAL (wp)           ::   zsec 
     59      REAL (wp)           ::   zfjulday 
    5660      !! 
    5761      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character 
    5862      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name 
    59       CHARACTER(lc)       ::   clpath   ! full path to ocean output restart file 
     63      CHARACTER(LEN=150)  ::   clpath   ! full path to ocean output restart file 
    6064      !!---------------------------------------------------------------------- 
    6165      ! 
    … …  
    8185      IF( kt == nitrst - 1 .OR. nstock == 1 .OR. ( kt == nitend .AND. .NOT. lrst_oce ) ) THEN 
    8286         IF( nitrst <= nitend .AND. nitrst > 0 ) THEN  
    83             ! beware of the format used to write kt (default is i8.8, that should be large enough...) 
    84             IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
    85             ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst 
     87            IF ( ln_rstdate ) THEN 
     88               zfjulday = fjulday + rdttra(1) / rday 
     89               IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday )   zfjulday = REAL(NINT(zfjulday),wp)   ! avoid truncation error 
     90               CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )            
     91               WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday 
     92            ELSE 
     93               ! beware of the format used to write kt (default is i8.8, that should be large enough...) 
     94               IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
     95               ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst 
     96               ENDIF 
    8697            ENDIF 
    8798            ! create the file 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r10246 r10249  
    121121   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sprecip           !: solid precipitation                          [Kg/m2/s] 
    122122   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fr_i              !: ice fraction = 1 - lead fraction      (between 0 to 1) 
     123   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   pressnow          !: UKMO SHELF pressure 
     124   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   apgu              !: UKMO SHELF pressure forcing 
     125   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   apgv              !: UKMO SHELF pressure forcing 
    123126#if defined key_cpl_carbon_cycle 
    124127   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   atm_co2           !: atmospheric pCO2                             [ppm] 
    … …  
    171174#endif 
    172175         &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) ,      & 
     176         &      pressnow(jpi,jpj), apgu(jpi,jpj)    , apgv(jpi,jpj) ,     & 
    173177         &      ssv_m  (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) ) 
    174178         ! 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90

    r10246 r10249  
    2222   USE lib_mpp         ! distribued memory computing library 
    2323   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     24   USE wrk_nemo        ! work arrays 
    2425 
    2526   IMPLICIT NONE 
    … …  
    2829   PUBLIC sbc_flx       ! routine called by step.F90 
    2930 
    30    INTEGER , PARAMETER ::   jpfld   = 5   ! maximum number of files to read  
     31   INTEGER , PARAMETER ::   jpfld   = 6   ! maximum number of files to read  
    3132   INTEGER , PARAMETER ::   jp_utau = 1   ! index of wind stress (i-component) file 
    3233   INTEGER , PARAMETER ::   jp_vtau = 2   ! index of wind stress (j-component) file 
    … …  
    3435   INTEGER , PARAMETER ::   jp_qsr  = 4   ! index of solar heat file 
    3536   INTEGER , PARAMETER ::   jp_emp  = 5   ! index of evaporation-precipation file 
     37   INTEGER , PARAMETER ::   jp_press = 6  ! index of pressure for UKMO shelf fluxes 
    3638   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read) 
     39   LOGICAL , PUBLIC    ::   ln_shelf_flx = .FALSE. ! UKMO SHELF specific flux flag 
     40   LOGICAL , PUBLIC    ::   ln_rel_wind = .FALSE.  ! UKMO SHELF specific flux flag - relative winds 
     41   REAL(wp)            ::   rn_wfac                ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem) 
     42   INTEGER             ::   jpfld_local   ! maximum number of files to read (locally modified depending on ln_shelf_flx)  
    3743 
    3844   !! * Substitutions 
    … …  
    8288      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient 
    8389      REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables 
     90      REAL     ::   cs                    ! UKMO SHELF: Friction co-efficient at surface 
     91      REAL     ::   totwindspd            ! UKMO SHELF: Magnitude of wind speed vector 
     92      REAL(wp), DIMENSION(:,:), POINTER ::   zwnd_i, zwnd_j    ! wind speed components at T-point 
     93 
     94      REAL(wp) ::   rhoa  = 1.22         ! Air density kg/m3 
     95      REAL(wp) ::   cdrag = 1.5e-3       ! drag coefficient  
    8496      !! 
    8597      CHARACTER(len=100) ::  cn_dir                               ! Root directory for location of flx files 
    8698      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                    ! array of namelist information structures 
    87       TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp  ! informations about the fields to be read 
    88       NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp 
     99      TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp, sn_press  !  informations about the fields to be read  
     100      LOGICAL     ::   ln_foam_flx  = .FALSE.                     ! UKMO FOAM specific flux flag  
     101      NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp,   &  
     102      &                    ln_foam_flx, sn_press, ln_shelf_flx, ln_rel_wind,    & 
     103      &                    rn_wfac 
    89104      !!--------------------------------------------------------------------- 
    90105      ! 
    … …  
    109124         slf_i(jp_emp ) = sn_emp 
    110125         ! 
    111          ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure 
     126            ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure 
     127            IF( ln_shelf_flx ) slf_i(jp_press) = sn_press 
     128    
     129            ! define local jpfld depending on shelf_flx logical 
     130            IF( ln_shelf_flx ) THEN 
     131               jpfld_local = jpfld 
     132            ELSE 
     133               jpfld_local = jpfld-1 
     134            ENDIF 
     135            ! 
    112136         IF( ierror > 0 ) THEN    
    113137            CALL ctl_stop( 'sbc_flx: unable to allocate sf structure' )   ;   RETURN   
    114138         ENDIF 
    115          DO ji= 1, jpfld 
     139         DO ji= 1, jpfld_local 
    116140            ALLOCATE( sf(ji)%fnow(jpi,jpj,1) ) 
    117141            IF( slf_i(ji)%ln_tint ) ALLOCATE( sf(ji)%fdta(jpi,jpj,1,2) ) 
    … …  
    128152      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN                        ! update ocean fluxes at each SBC frequency 
    129153 
     154         !!UKMO SHELF wind speed relative to surface currents - put here to allow merging with coupling branch 
     155         IF( ln_shelf_flx ) THEN 
     156            CALL wrk_alloc( jpi,jpj, zwnd_i, zwnd_j ) 
     157 
     158            IF( ln_rel_wind ) THEN 
     159               DO jj = 1, jpj 
     160                  DO ji = 1, jpi 
     161                     zwnd_i(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) - rn_wfac * ssu_m(ji,jj) 
     162                     zwnd_j(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) - rn_wfac * ssv_m(ji,jj) 
     163                  END DO 
     164               END DO 
     165            ELSE 
     166               zwnd_i(:,:) = sf(jp_utau)%fnow(:,:,1) 
     167               zwnd_j(:,:) = sf(jp_vtau)%fnow(:,:,1) 
     168            ENDIF 
     169         ENDIF 
     170 
    130171         IF( ln_dm2dc ) THEN   ;   qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )   ! modify now Qsr to include the diurnal cycle 
    131172         ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1) 
    132173         ENDIF 
    133174!CDIR COLLAPSE 
     175            !!UKMO SHELF effect of atmospheric pressure on SSH 
     176            ! If using ln_apr_dyn, this is done there so don't repeat here. 
     177            IF( ln_shelf_flx .AND. .NOT. ln_apr_dyn) THEN 
     178               DO jj = 1, jpjm1 
     179                  DO ji = 1, jpim1 
     180                     apgu(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji+1,jj,1)-sf(jp_press)%fnow(ji,jj,1))/e1u(ji,jj) 
     181                     apgv(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji,jj+1,1)-sf(jp_press)%fnow(ji,jj,1))/e2v(ji,jj) 
     182                  END DO 
     183               END DO 
     184            ENDIF ! ln_shelf_flx 
     185 
    134186         DO jj = 1, jpj                                           ! set the ocean fluxes from read fields 
    135187            DO ji = 1, jpi 
    136                utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) 
    137                vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) 
    138                qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) 
    139                emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1) 
     188                   IF( ln_shelf_flx ) THEN 
     189                      !! UKMO SHELF - need atmospheric pressure to calculate Haney forcing 
     190                      pressnow(ji,jj) = sf(jp_press)%fnow(ji,jj,1) 
     191                      !! UKMO SHELF flux files contain wind speed not wind stress 
     192                      totwindspd = sqrt(zwnd_i(ji,jj)*zwnd_i(ji,jj) + zwnd_j(ji,jj)*zwnd_j(ji,jj)) 
     193                      cs = 0.63 + (0.066 * totwindspd) 
     194                      utau(ji,jj) = cs * (rhoa/rau0) * zwnd_i(ji,jj) * totwindspd 
     195                      vtau(ji,jj) = cs * (rhoa/rau0) * zwnd_j(ji,jj) * totwindspd 
     196                   ELSE 
     197                      utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) 
     198                      vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) 
     199                   ENDIF 
     200                   qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1) 
     201                   IF( ln_foam_flx .OR. ln_shelf_flx ) THEN 
     202                      !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot) 
     203                      qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) 
     204                      !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P 
     205                      emp (ji,jj) = -1. * sf(jp_emp )%fnow(ji,jj,1) 
     206                   ELSE 
     207                      qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) 
     208                      emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1) 
     209                   ENDIF 
    140210            END DO 
    141211         END DO 
    … …  
    143213         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST 
    144214         ! 
     215    
     216            !! UKMO FOAM wind fluxes need lbc_lnk calls owing to a bug in interp.exe 
     217            IF( ln_foam_flx ) THEN 
     218               CALL lbc_lnk( utau(:,:), 'U', -1. ) 
     219               CALL lbc_lnk( vtau(:,:), 'V', -1. ) 
     220            ENDIF 
     221     
    145222         !                                                        ! module of wind stress and wind speed at T-point 
    146223         zcoef = 1. / ( zrhoa * zcdrag ) 
    … …  
    162239            WRITE(numout,*)  
    163240            WRITE(numout,*) '        read daily momentum, heat and freshwater fluxes OK' 
    164             DO jf = 1, jpfld 
     241            DO jf = 1, jpfld_local 
    165242               IF( jf == jp_utau .OR. jf == jp_vtau )   zfact =     1. 
    166243               IF( jf == jp_qtot .OR. jf == jp_qsr  )   zfact =     0.1 
    … …  
    173250         ENDIF 
    174251         ! 
     252         IF( ln_shelf_flx ) THEN 
     253            CALL wrk_dealloc( jpi,jpj, zwnd_i, zwnd_j ) 
     254         ENDIF 
     255         ! 
    175256      ENDIF 
    176257      ! 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90

    r10246 r10249  
    4242   LOGICAL         ::   ln_sssr_bnd     ! flag to bound erp term  
    4343   REAL(wp)        ::   rn_sssr_bnd     ! ABS(Max./Min.) value of erp term [mm/day] 
     44   LOGICAL         ::   ln_UKMO_haney   ! UKMO specific flag to calculate Haney forcing   
    4445 
    4546   REAL(wp) , ALLOCATABLE, DIMENSION(:) ::   buffer   ! Temporary buffer for exchange 
    … …  
    7980      INTEGER  ::   ierror   ! return error code 
    8081      !! 
     82      REAL(wp) ::   sst1,sst2                      ! sea surface temperature 
     83      REAL(wp) ::   e_sst1, e_sst2                 ! saturation vapour pressure 
     84      REAL(wp) ::   qs1,qs2                        ! specific humidity 
     85      REAL(wp) ::   pr_tmp                         ! temporary variable for pressure 
     86  
     87      REAL(wp), DIMENSION(jpi,jpj) ::  hny_frc1    ! Haney forcing for sensible heat, correction for latent heat    
     88      REAL(wp), DIMENSION(jpi,jpj) ::  hny_frc2    ! Haney forcing for sensible heat, correction for latent heat    
     89      !! 
    8190      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files 
    8291      TYPE(FLD_N) ::   sn_sst, sn_sss        ! informations about the fields to be read 
    … …  
    95104            ! 
    96105            IF( nn_sstr == 1 ) THEN                                   !* Temperature restoring term 
    97                DO jj = 1, jpj 
    98                   DO ji = 1, jpi 
    99                      zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) ) 
    100                      qns(ji,jj) = qns(ji,jj) + zqrp 
    101                      qrp(ji,jj) = zqrp 
    102                   END DO 
    103                END DO 
     106                  IF( ln_UKMO_haney ) THEN 
     107                     DO jj = 1, jpj 
     108                        DO ji = 1, jpi 
     109                           sst1   =  sst_m(ji,jj) 
     110                           sst2   =  sf_sst(1)%fnow(ji,jj,1)    
     111                           e_sst1 = 10**((0.7859+0.03477*sst1)/(1.+0.00412*sst1)) 
     112                           e_sst2 = 10**((0.7859+0.03477*sst2)/(1.+0.00412*sst2))          
     113                           pr_tmp = 0.01*pressnow(ji,jj)  !pr_tmp = 1012.0 
     114                           qs1    = (0.62197*e_sst1)/(pr_tmp-0.378*e_sst1) 
     115                           qs2    = (0.62197*e_sst2)/(pr_tmp-0.378*e_sst2) 
     116                           hny_frc1(ji,jj) = sst1-sst2                    
     117                           hny_frc2(ji,jj) = qs1-qs2                      
     118                          !Might need to mask off land points. 
     119                           hny_frc1(ji,jj)=-hny_frc1(ji,jj)*wndm(ji,jj)*1.42 
     120                           hny_frc2(ji,jj)=-hny_frc2(ji,jj)*wndm(ji,jj)*4688.0 
     121                           qns(ji,jj)=qns(ji,jj)+hny_frc1(ji,jj)+hny_frc2(ji,jj)    
     122                           qrp(ji,jj) = 0.e0 
     123                        END DO 
     124                     END DO 
     125                  ELSE 
     126                     DO jj = 1, jpj 
     127                        DO ji = 1, jpi 
     128                           zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) ) 
     129                           qns(ji,jj) = qns(ji,jj) + zqrp 
     130                           qrp(ji,jj) = zqrp 
     131                        END DO 
     132                     END DO 
     133                  ENDIF 
    104134               CALL iom_put( "qrp", qrp )                             ! heat flux damping 
    105135            ENDIF 
    … …  
    163193      CHARACTER(len=100) ::  cn_dir          ! Root directory for location of ssr files 
    164194      TYPE(FLD_N) ::   sn_sst, sn_sss        ! informations about the fields to be read 
    165       NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, sn_sss, ln_sssr_bnd, rn_sssr_bnd 
     195      NAMELIST/namsbc_ssr/ cn_dir, nn_sstr, nn_sssr, rn_dqdt, rn_deds, sn_sst, sn_sss, ln_sssr_bnd, rn_sssr_bnd, ln_UKMO_haney 
    166196      INTEGER     ::  ios 
    167197      !!---------------------------------------------------------------------- 
    … …  
    189219         WRITE(numout,*) '      flag to bound erp term                 ln_sssr_bnd = ', ln_sssr_bnd 
    190220         WRITE(numout,*) '      ABS(Max./Min.) erp threshold           rn_sssr_bnd = ', rn_sssr_bnd, ' mm/day' 
     221         WRITE(numout,*) '      Haney forcing                          ln_UKMO_haney = ', ln_UKMO_haney 
    191222      ENDIF 
    192223      ! 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/SBC/tide.h90

    r4292 r10249  
    2828   Wave(18) = tide(  'L2'     , 0.006694 ,    2   ,  2 , -1 ,  2 , -1 ,  0  , +180  ,  2   , -2   ,  0   ,  0   , 0 ,  215    ) 
    2929   Wave(19) = tide(  'T2'     , 0.006614 ,    2   ,  2 ,  0 , -1 ,  0 ,  1  ,    0  ,  0   ,  0   ,  0   ,  0   , 0 ,    0    ) 
     30   ! 
     31   !             !! name_tide , equitide , nutide , nt , ns , nh , np , np1 , shift , nksi , nnu0 , nnu1 , nnu2 , R , formula !! 
     32   Wave(20) = tide(  'MNS2'   , 0.000000 ,    2   ,  2 , -5 ,  4 ,  1 ,  0  ,    0  ,  4   , -4   ,  0   ,  0   , 0 ,    6    ) 
     33   Wave(21) = tide(  'Lam2'   , 0.001760 ,    2   ,  2 , -1 ,  0 ,  1 ,  0  , +180  ,  2   , -2   ,  0   ,  0   , 0 ,   78    ) 
     34   Wave(22) = tide(  'MSN2'   , 0.000000 ,    2   ,  2 ,  1 ,  0 ,  1 ,  0  ,    0  ,  2   , -2   ,  0   ,  2   , 0 ,    6    ) 
     35   Wave(23) = tide(  '2SM2'   , 0.000000 ,    2   ,  2 ,  2 , -2 ,  0 ,  0  ,    0  , -2   ,  2   ,  0   ,  0   , 0 ,   16    ) 
     36   Wave(24) = tide(  'MO3'    , 0.000000 ,    3   ,  3 , -4 ,  1 ,  0 ,  0  ,  +90  ,  2   , -2   ,  0   ,  0   , 0 ,   13    ) 
     37   Wave(25) = tide(  'MK3'    , 0.000000 ,    3   ,  3 , -2 ,  3 ,  0 ,  0  ,  -90  ,  2   , -2   , -1   ,  0   , 0 ,   10    ) 
     38   Wave(26) = tide(  'MN4'    , 0.000000 ,    4   ,  4 , -5 ,  4 ,  1 ,  0  ,    0  ,  4   , -4   ,  0   ,  0   , 0 ,    1    ) 
     39   Wave(27) = tide(  'MS4'    , 0.000000 ,    4   ,  4 , -2 ,  2 ,  0 ,  0  ,    0  ,  2   , -2   ,  0   ,  0   , 0 ,    2    ) 
     40   Wave(28) = tide(  'M6'     , 0.000000 ,    6   ,  6 , -6 ,  6 ,  0 ,  0  ,    0  ,  6   , -6   ,  0   ,  0   , 0 ,    4    ) 
     41   Wave(29) = tide(  '2MS6'   , 0.000000 ,    6   ,  6 , -4 ,  4 ,  0 ,  0  ,    0  ,  4   , -4   ,  0   ,  0   , 0 ,    6    ) 
     42   Wave(30) = tide(  '2MK6'   , 0.000000 ,    6   ,  6 , -4 ,  6 ,  0 ,  0  ,    0  ,  4   , -4   ,  0   , -2   , 0 ,    5    ) 
     43   Wave(31) = tide(  '3M2S2'  , 0.000000 ,    2   , 2  , -6 ,  6 ,  0 ,  0  ,    0  ,  6   , -6   ,  0   ,  0   , 0 ,   12    ) 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/SBC/tide_mod.F90

    r10246 r10249  
    1616   PUBLIC   tide_init_Wave   ! called by tideini and diaharm modules 
    1717 
    18    INTEGER, PUBLIC, PARAMETER ::   jpmax_harmo = 19   !: maximum number of harmonic 
     18   INTEGER, PUBLIC, PARAMETER ::   jpmax_harmo = 31   !: maximum number of harmonic 
    1919 
    2020   TYPE, PUBLIC ::    tide 
    21       CHARACTER(LEN=4) ::   cname_tide 
     21      CHARACTER(LEN=5) ::   cname_tide 
    2222      REAL(wp)         ::   equitide 
    2323      INTEGER          ::   nutide 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r10248 r10249  
    12411241      IF(lwm) WRITE( numond, nameos ) 
    12421242      ! 
    1243       rau0        = 1026._wp                 !: volumic mass of reference     [kg/m3] 
     1243      rau0        = 1020._wp                 !: volumic mass of reference     [kg/m3] 
     1244!     rau0        = 1026._wp                 !: volumic mass of reference     [kg/m3] 
    12441245      rcp         = 3991.86795711963_wp      !: heat capacity     [J/K] 
    12451246      ! 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90

    r10248 r10249  
    100100         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 
    101101      ENDIF 
     102! slwa unless you use l_trdtra too, the above switches off trend calculations for l_trdtrc 
     103         l_trd = .FALSE. 
     104         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 
     105!slwa 
    102106      ! 
    103107      IF( l_trd )  THEN 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

    r10248 r10249  
    5858 
    5959   REAL(wp) ::   rbcp   ! Brown & Campana parameters for semi-implicit hpg 
     60   INTEGER  ::   warn_1, warn_2   ! indicators for warning statement 
    6061 
    6162   !! * Substitutions 
    … …  
    9394      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    9495      !! 
    95       INTEGER  ::   jk, jn    ! dummy loop indices 
    96       REAL(wp) ::   zfact     ! local scalars 
     96      INTEGER  ::   jk, jn, ji, jj     ! dummy loop indices 
     97      REAL(wp) ::   zfact, zfreeze     ! local scalars 
    9798      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds 
    9899      !!---------------------------------------------------------------------- 
    … …  
    125126      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dtra(:) = 2._wp* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog) 
    126127      ENDIF 
     128 
     129#if ( ! defined key_lim3 && ! defined key_lim2 && ! key_cice ) 
     130      IF ( kt == nit000 ) warn_1=0 
     131      warn_2=0 
     132      DO jk = 1, jpkm1 
     133         DO jj = 1, jpj 
     134            DO ji = 1, jpi 
     135               IF ( tsa(ji,jj,jk,jp_tem) .lt. 0.0 ) THEN 
     136                  ! calculate freezing point 
     137                  zfreeze = ( -0.0575_wp + 1.710523E-3 * Sqrt(Abs(tsn(ji,jj,jk,jp_sal)))   &  
     138                            - 2.154996E-4 * tsn(ji,jj,jk,jp_sal) ) * tsn(ji,jj,jk,jp_sal) - 7.53E-4 * ( 10.0_wp + fsdept(ji,jj,jk) ) 
     139                  IF ( tsa(ji,jj,jk,jp_tem) .lt. zfreeze ) THEN 
     140                     tsa(ji,jj,jk,jp_tem)=zfreeze 
     141                     warn_2=1 
     142                  ENDIF 
     143               ENDIF 
     144            END DO 
     145         END DO 
     146      END DO 
     147      CALL mpp_max(warn_1) 
     148      CALL mpp_max(warn_2) 
     149      IF ( (warn_1 == 0) .and. (warn_2 /= 0) ) THEN 
     150         IF(lwp) THEN 
     151            CALL ctl_warn( ' Temperatures dropping below freezing point, ', & 
     152                      &    ' being forced to freezing point, no longer conservative' )  
     153         ENDIF 
     154         warn_1=1 
     155      ENDIF 
     156#endif 
    127157 
    128158      ! trends computation initialisation 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r10246 r10249  
    4646   LOGICAL , PUBLIC ::   ln_qsr_ice   !: light penetration for ice-model LIM3 (clem) 
    4747   INTEGER , PUBLIC ::   nn_chldta    !: use Chlorophyll data (=1) or not (=0) 
     48   INTEGER , PUBLIC ::   nn_kd490dta  !: use kd490dta data (=1) or not (=0) 
    4849   REAL(wp), PUBLIC ::   rn_abs       !: fraction absorbed in the very near surface (RGB & 2 bands) 
    4950   REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands) 
    … …  
    5455   REAL(wp) ::   xsi1r                           !: inverse of rn_si1 
    5556   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read) 
     57   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_kd490 ! structure of input kd490 (file informations, fields read) 
    5658   INTEGER, PUBLIC ::   nksr              ! levels below which the light cannot penetrate ( depth larger than 391 m) 
    5759   REAL(wp), DIMENSION(3,61) ::   rkrgb   !: tabulated attenuation coefficients for RGB absorption 
    … …  
    306308            ! 
    307309         ENDIF 
     310! slwa 
     311         IF( nn_kd490dta == 1 ) THEN                      !  use KD490 data read in   ! 
     312            !                                             ! ------------------------- ! 
     313               nksr = jpk - 1 
     314               ! 
     315               CALL fld_read( kt, 1, sf_kd490 )     ! Read kd490 data and provide it at the current time step 
     316               ! 
     317               zcoef  = ( 1. - rn_abs ) 
     318               ze0(:,:,1) = rn_abs  * qsr(:,:) 
     319               ze1(:,:,1) = zcoef * qsr(:,:) 
     320               zea(:,:,1) =         qsr(:,:) 
     321               ! 
     322               DO jk = 2, nksr+1 
     323!CDIR NOVERRCHK 
     324                  DO jj = 1, jpj 
     325!CDIR NOVERRCHK    
     326                     DO ji = 1, jpi 
     327                        zc0 = ze0(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * xsi0r     ) 
     328                        zc1 = ze1(ji,jj,jk-1) * EXP( - fse3t(ji,jj,jk-1) * sf_kd490(1)%fnow(ji,jj,1) ) 
     329                        ze0(ji,jj,jk) = zc0 
     330                        ze1(ji,jj,jk) = zc1 
     331                        zea(ji,jj,jk) = ( zc0 + zc1 ) * tmask(ji,jj,jk) 
     332                     END DO 
     333                  END DO 
     334               END DO 
     335               ! clem: store attenuation coefficient of the first ocean level 
     336               IF ( ln_qsr_ice ) THEN 
     337                  DO jj = 1, jpj 
     338                     DO ji = 1, jpi 
     339                        zzc0 = rn_abs * EXP( - fse3t(ji,jj,1) * xsi0r     ) 
     340                        zzc1 = zcoef  * EXP( - fse3t(ji,jj,1) * sf_kd490(1)%fnow(ji,jj,1) ) 
     341                        fraqsr_1lev(ji,jj) = 1.0 - ( zzc0 + zzc1 ) * tmask(ji,jj,2)  
     342                     END DO 
     343                  END DO 
     344               ENDIF 
     345               ! 
     346               DO jk = 1, nksr                                        ! compute and add qsr trend to ta 
     347                  qsr_hc(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) ) 
     348               END DO 
     349               zea(:,:,nksr+1:jpk) = 0.e0     !  
     350               CALL iom_put( 'qsr3d', zea )   ! Shortwave Radiation 3D distribution 
     351               ! 
     352        ENDIF   ! use KD490 data 
     353!slwa 
    308354         ! 
    309355         !                                        Add to the general trend 
    … …  
    374420      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files 
    375421      TYPE(FLD_N)        ::   sn_chl   ! informations about the chlorofyl field to be read 
    376       !! 
    377       NAMELIST/namtra_qsr/  sn_chl, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  & 
    378          &                  nn_chldta, rn_abs, rn_si0, rn_si1 
     422      TYPE(FLD_N)        ::   sn_kd490 ! informations about the kd490 field to be read 
     423      !! 
     424      NAMELIST/namtra_qsr/  sn_chl, sn_kd490, cn_dir, ln_traqsr, ln_qsr_rgb, ln_qsr_2bd, ln_qsr_bio, ln_qsr_ice,  & 
     425         &                  nn_chldta, rn_abs, rn_si0, rn_si1, nn_kd490dta 
    379426      !!---------------------------------------------------------------------- 
    380427 
    … …  
    409456         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0 = ', rn_si0 
    410457         WRITE(numout,*) '      2 bands: longest depth of extinction         rn_si1 = ', rn_si1 
     458         WRITE(numout,*) '      read in KD490 data                       nn_kd490dta  = ', nn_kd490dta 
    411459      ENDIF 
    412460 
    … …  
    422470         IF( ln_qsr_2bd  )   ioptio = ioptio + 1 
    423471         IF( ln_qsr_bio  )   ioptio = ioptio + 1 
     472         IF( nn_kd490dta == 1 )   ioptio = ioptio + 1 
    424473         ! 
    425474         IF( ioptio /= 1 ) & 
    … …  
    431480         IF( ln_qsr_2bd                      )   nqsr =  3 
    432481         IF( ln_qsr_bio                      )   nqsr =  4 
     482         IF( nn_kd490dta == 1                )   nqsr =  5 
    433483         ! 
    434484         IF(lwp) THEN                   ! Print the choice 
    … …  
    438488            IF( nqsr ==  3 )   WRITE(numout,*) '         2 bands light penetration' 
    439489            IF( nqsr ==  4 )   WRITE(numout,*) '         bio-model light penetration' 
     490            IF( nqsr ==  5 )   WRITE(numout,*) '         KD490 light penetration' 
    440491         ENDIF 
    441492         ! 
    … …  
    447498         xsi0r = 1.e0 / rn_si0 
    448499         xsi1r = 1.e0 / rn_si1 
     500         IF( nn_kd490dta == 1 ) THEN           !* KD490 data : set sf_kd490 structure 
     501            IF(lwp) WRITE(numout,*) 
     502            IF(lwp) WRITE(numout,*) '        KD490 read in a file' 
     503            ALLOCATE( sf_kd490(1), STAT=ierror ) 
     504            IF( ierror > 0 ) THEN 
     505               CALL ctl_stop( 'tra_qsr_init: unable to allocate sf_kd490 structure' )   ;   RETURN 
     506            ENDIF 
     507            ALLOCATE( sf_kd490(1)%fnow(jpi,jpj,1)   ) 
     508            IF( sn_kd490%ln_tint )ALLOCATE( sf_kd490(1)%fdta(jpi,jpj,1,2) ) 
     509            !                                        ! fill sf_kd490 with sn_kd490 and control print 
     510            CALL fld_fill( sf_kd490, (/ sn_kd490 /), cn_dir, 'tra_qsr_init',   & 
     511               &                                         'Solar penetration function of read KD490', 'namtra_qsr' ) 
    449512         !                                ! ---------------------------------- ! 
    450          IF( ln_qsr_rgb ) THEN            !  Red-Green-Blue light penetration  ! 
     513         ELSEIF( ln_qsr_rgb ) THEN            !  Red-Green-Blue light penetration  ! 
    451514            !                             ! ---------------------------------- ! 
    452515            ! 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r10248 r10249  
    2525   USE trd_oce         ! trends: ocean variables 
    2626   USE trdtra          ! trends manager: tracers  
     27   USE tradwl          ! solar radiation penetration (downwell method) 
    2728   ! 
    2829   USE in_out_manager  ! I/O manager 
    … …  
    3334   USE timing          ! Timing 
    3435   USE eosbn2 
     36#if defined key_asminc    
     37   USE asminc          ! Assimilation increment 
     38#endif 
    3539 
    3640   IMPLICIT NONE 
    … …  
    138142 
    139143!!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration 
    140       IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration 
     144      IF( .NOT.ln_traqsr .and. .NOT.ln_tradwl ) THEN     ! no solar radiation penetration 
    141145         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns 
    142146         qsr(:,:) = 0.e0                     ! qsr set to zero 
    … …  
    278282         END DO   
    279283      ENDIF 
     284 
     285#if defined key_asminc 
     286! WARNING: THIS MAY WELL NOT BE REQUIRED - WE DON'T WANT TO CHANGE T&S BUT THIS MAY COMPENSATE ANOTHER TERM... 
     287! Rate of change in e3t for each level is ssh_iau*e3t_0/ht_0 
     288! Contribution to tsa should be rate of change in level / per m of ocean? (hence the division by fse3t_n) 
     289      IF( ln_sshinc ) THEN         ! input of heat and salt due to assimilation 
     290         DO jj = 2, jpj  
     291            DO ji = fs_2, fs_jpim1 
     292               zdep = ssh_iau(ji,jj) / ( ht_0(ji,jj) + 1.0 - ssmask(ji, jj) ) 
     293               DO jk = 1, jpkm1 
     294                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)   & 
     295                                        &            + tsn(ji,jj,jk,jp_tem) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) ) 
     296                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)   & 
     297                                        &            + tsn(ji,jj,jk,jp_sal) * zdep * ( e3t_0(ji,jj,jk) / fse3t_n(ji,jj,jk) ) 
     298               END DO 
     299            END DO   
     300         END DO   
     301      ENDIF 
     302#endif 
    280303  
    281304      IF( l_trdtra )   THEN                      ! send trends for further diagnostics 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90

    r10246 r10249  
    203203         DO jj = 2, jpjm1 
    204204            DO ji = fs_2, fs_jpim1   ! vector opt. 
     205#if defined key_tracer_budget 
     206!              ptrd(ji,jj,jk) = - (     pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik)  )  * tmask(ji,jj,jk) 
     207               ptrd(ji,jj,jk) = -      pf (ji,jj,jk) * tmask(ji,jj,jk) 
     208#else 
    205209               ptrd(ji,jj,jk) = - (     pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik)                        & 
    206210                 &                  - ( pun(ji,jj,jk) - pun(ji-ii,jj-ij,jk-ik) ) * ptn(ji,jj,jk)  )   & 
    207211                 &              / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )  * tmask(ji,jj,jk) 
     212#endif 
    208213            END DO 
    209214         END DO 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

    r10248 r10249  
    1818   USE phycst          ! physical constants 
    1919   USE iom             ! I/O library 
     20   USE eosbn2          ! for zdf_mxl_zint 
    2021   USE lib_mpp         ! MPP library 
    2122   USE wrk_nemo        ! work arrays 
    … …  
    2627   PRIVATE 
    2728 
     29   PUBLIC   zdf_mxl_tref  ! called by asminc.F90 
    2830   PUBLIC   zdf_mxl       ! called by step.F90 
    29    PUBLIC   zdf_mxl_alloc ! Used in zdf_tke_init 
    30  
     31 
     32   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld_tref  !: mixed layer depth at t-points - temperature criterion [m] 
    3133   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln    !: number of level in the mixed layer (used by TOP) 
    3234   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld    !: mixing layer depth (turbocline)      [m] 
    3335   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m] 
    3436   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m] 
     37   REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:) ::   hmld_zint  !: vertically-interpolated mixed layer depth   [m]  
     38   REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:) ::   htc_mld    ! Heat content of hmld_zint 
     39   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    :: ll_found   ! Is T_b to be found by interpolation ?  
     40   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: ll_belowml ! Flag points below mixed layer when ll_found=F 
    3541 
    3642   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
    3743   REAL(wp)         ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth 
     44 
     45   TYPE, PUBLIC :: MXL_ZINT   !: Structure for MLD defs 
     46      INTEGER   :: mld_type   ! mixed layer type      
     47      REAL(wp)  :: zref       ! depth of initial T_ref 
     48      REAL(wp)  :: dT_crit    ! Critical temp diff 
     49      REAL(wp)  :: iso_frac   ! Fraction of rn_dT_crit used 
     50   END TYPE MXL_ZINT 
     51 
     52!Used for 25h mean 
     53   LOGICAL, PRIVATE :: mld_25h_init = .TRUE.    !Logical used to initalise 25h 
     54                                                !outputs. Necassary, because we need to 
     55                                                !initalise the mld_25h on the zeroth 
     56                                                !timestep (i.e in the nemogcm_init call) 
     57   LOGICAL, PRIVATE :: mld_25h_write = .FALSE.  !Logical confirm 25h calculating/processing 
     58 
     59   INTEGER, SAVE :: i_cnt_25h                   ! Counter for 25 hour means 
     60 
     61   REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   hmld_zint_25h 
    3862 
    3963   !! * Substitutions 
    … …  
    5276      zdf_mxl_alloc = 0      ! set to zero if no array to be allocated 
    5377      IF( .NOT. ALLOCATED( nmln ) ) THEN 
    54          ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc ) 
     78         ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj),       & 
     79        &          htc_mld(jpi,jpj),                                                                    & 
     80        &          ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc ) 
    5581         ! 
     82         ALLOCATE(hmld_tref(jpi,jpj)) 
    5683         IF( lk_mpp             )   CALL mpp_sum ( zdf_mxl_alloc ) 
    5784         IF( zdf_mxl_alloc /= 0 )   CALL ctl_warn('zdf_mxl_alloc: failed to allocate arrays.') 
    … …  
    6087   END FUNCTION zdf_mxl_alloc 
    6188 
     89   SUBROUTINE zdf_mxl_tref()  
     90      !!----------------------------------------------------------------------  
     91      !!                  ***  ROUTINE zdf_mxl_tref  ***  
     92      !!                     
     93      !! ** Purpose :   Compute the mixed layer depth with temperature criteria.  
     94      !!  
     95      !! ** Method  :   The temperature-defined mixed layer depth is required  
     96      !!                   when assimilating SST in a 2D analysis.   
     97      !!  
     98      !! ** Action  :   hmld_tref  
     99      !!----------------------------------------------------------------------  
     100      !  
     101      INTEGER  ::   ji, jj, jk   ! dummy loop indices  
     102      REAL(wp) ::   t_ref               ! Reference temperature    
     103      REAL(wp) ::   temp_c = 0.2        ! temperature criterion for mixed layer depth    
     104      !!----------------------------------------------------------------------  
     105      !  
     106      ! Initialise array  
     107      IF( zdf_mxl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl_tref : unable to allocate arrays' )  
     108        
     109      !For the AMM model assimiation uses a temperature based mixed layer depth    
     110      !This is defined here    
     111      DO jj = 1, jpj    
     112         DO ji = 1, jpi    
     113           hmld_tref(ji,jj)=fsdept(ji,jj,1  )     
     114           IF(ssmask(ji,jj) > 0.)THEN    
     115             t_ref=tsn(ji,jj,1,jp_tem)   
     116             DO jk=2,jpk    
     117               IF(ssmask(ji,jj)==0.)THEN    
     118                  hmld_tref(ji,jj)=fsdept(ji,jj,jk )    
     119                  EXIT    
     120               ELSEIF( ABS(tsn(ji,jj,jk,jp_tem)-t_ref) < temp_c)THEN    
     121                  hmld_tref(ji,jj)=fsdept(ji,jj,jk )    
     122               ELSE    
     123                  EXIT    
     124               ENDIF    
     125             ENDDO    
     126           ENDIF    
     127         ENDDO    
     128      ENDDO  
     129    
     130   END SUBROUTINE zdf_mxl_tref  
    62131 
    63132   SUBROUTINE zdf_mxl( kt ) 
    … …  
    128197            iikn = nmln(ji,jj) 
    129198            imkt = mikt(ji,jj) 
    130             hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj)    ! Turbocline depth  
    131             hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj)    ! Mixed layer depth 
    132             hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer 
     199            hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,imkt )            )  * ssmask(ji,jj)    ! Turbocline depth  
     200            hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,MAX( imkt,nla10 ) ) ) * ssmask(ji,jj)    ! Mixed layer depth 
     201            hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt )            )  * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer 
    133202         END DO 
    134203      END DO 
    … …  
    138207      ENDIF 
    139208       
     209      ! Vertically-interpolated mixed-layer depth diagnostic 
     210      CALL zdf_mxl_zint( kt ) 
     211 
    140212      IF(ln_ctl)   CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 ) 
    141213      ! 
    … …  
    145217      ! 
    146218   END SUBROUTINE zdf_mxl 
     219 
     220   SUBROUTINE zdf_mxl_zint_mld( sf )  
     221      !!----------------------------------------------------------------------------------  
     222      !!                    ***  ROUTINE zdf_mxl_zint_mld  ***  
     223      !                                                                         
     224      !   Calculate vertically-interpolated mixed layer depth diagnostic.  
     225      !             
     226      !   This routine can calculate the mixed layer depth diagnostic suggested by 
     227      !   Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate 
     228      !   vertically-interpolated mixed-layer depth diagnostics with other parameter 
     229      !   settings set in the namzdf_mldzint namelist.   
     230      !  
     231      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the   
     232      !   density has increased by an amount equivalent to a temperature difference of   
     233      !   0.8C at the surface.  
     234      !  
     235      !   For other values of mld_type the mixed layer is calculated as the depth at   
     236      !   which the temperature differs by 0.8C from the surface temperature.   
     237      !                                                                         
     238      !   David Acreman, Daley Calvert                                       
     239      !  
     240      !!-----------------------------------------------------------------------------------  
     241 
     242      TYPE(MXL_ZINT), INTENT(in)  :: sf 
     243 
     244      ! Diagnostic criteria 
     245      INTEGER   :: nn_mld_type   ! mixed layer type      
     246      REAL(wp)  :: rn_zref       ! depth of initial T_ref 
     247      REAL(wp)  :: rn_dT_crit    ! Critical temp diff 
     248      REAL(wp)  :: rn_iso_frac   ! Fraction of rn_dT_crit used 
     249 
     250      ! Local variables 
     251      REAL(wp), PARAMETER :: zepsilon = 1.e-30          ! local small value 
     252      INTEGER, POINTER, DIMENSION(:,:) :: ikmt          ! number of active tracer levels  
     253      INTEGER, POINTER, DIMENSION(:,:) :: ik_ref        ! index of reference level  
     254      INTEGER, POINTER, DIMENSION(:,:) :: ik_iso        ! index of last uniform temp level  
     255      REAL, POINTER, DIMENSION(:,:,:)  :: zT            ! Temperature or density  
     256      REAL, POINTER, DIMENSION(:,:)    :: ppzdep        ! depth for use in calculating d(rho)  
     257      REAL, POINTER, DIMENSION(:,:)    :: zT_ref        ! reference temperature  
     258      REAL    :: zT_b                                   ! base temperature  
     259      REAL, POINTER, DIMENSION(:,:,:)  :: zdTdz         ! gradient of zT  
     260      REAL, POINTER, DIMENSION(:,:,:)  :: zmoddT        ! Absolute temperature difference  
     261      REAL    :: zdz                                    ! depth difference  
     262      REAL    :: zdT                                    ! temperature difference  
     263      REAL, POINTER, DIMENSION(:,:)    :: zdelta_T      ! difference critereon  
     264      REAL, POINTER, DIMENSION(:,:)    :: zRHO1, zRHO2  ! Densities  
     265      INTEGER :: ji, jj, jk                             ! loop counter  
     266 
     267      !!-------------------------------------------------------------------------------------  
     268      !   
     269      CALL wrk_alloc( jpi, jpj, ikmt, ik_ref, ik_iso)  
     270      CALL wrk_alloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 )  
     271      CALL wrk_alloc( jpi, jpj, jpk, zT, zdTdz, zmoddT )  
     272 
     273      ! Unpack structure 
     274      nn_mld_type = sf%mld_type 
     275      rn_zref     = sf%zref 
     276      rn_dT_crit  = sf%dT_crit 
     277      rn_iso_frac = sf%iso_frac 
     278 
     279      ! Set the mixed layer depth criterion at each grid point  
     280      IF( nn_mld_type == 0 ) THEN 
     281         zdelta_T(:,:) = rn_dT_crit 
     282         zT(:,:,:) = rhop(:,:,:) 
     283      ELSE IF( nn_mld_type == 1 ) THEN 
     284         ppzdep(:,:)=0.0  
     285         call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) )  
     286! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST  
     287! [assumes number of tracers less than number of vertical levels]  
     288         zT(:,:,1:jpts)=tsn(:,:,1,1:jpts)  
     289         zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit  
     290         CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) )  
     291         zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0  
     292         ! RHO from eos (2d version) doesn't calculate north or east halo:  
     293         CALL lbc_lnk( zdelta_T, 'T', 1. )  
     294         zT(:,:,:) = rhop(:,:,:)  
     295      ELSE  
     296         zdelta_T(:,:) = rn_dT_crit                       
     297         zT(:,:,:) = tsn(:,:,:,jp_tem)                            
     298      END IF  
     299 
     300      ! Calculate the gradient of zT and absolute difference for use later  
     301      DO jk = 1 ,jpk-2  
     302         zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / fse3w(:,:,jk+1)  
     303         zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )  
     304      END DO  
     305 
     306      ! Find density/temperature at the reference level (Kara et al use 10m).           
     307      ! ik_ref is the index of the box centre immediately above or at the reference level  
     308      ! Find rn_zref in the array of model level depths and find the ref     
     309      ! density/temperature by linear interpolation.                                    
     310      DO jk = jpkm1, 2, -1  
     311         WHERE ( fsdept(:,:,jk) > rn_zref )  
     312           ik_ref(:,:) = jk - 1  
     313           zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - fsdept(:,:,jk-1) )  
     314         END WHERE  
     315      END DO  
     316 
     317      ! If the first grid box centre is below the reference level then use the  
     318      ! top model level to get zT_ref  
     319      WHERE ( fsdept(:,:,1) > rn_zref )   
     320         zT_ref = zT(:,:,1)  
     321         ik_ref = 1  
     322      END WHERE  
     323 
     324      ! The number of active tracer levels is 1 less than the number of active w levels  
     325      ikmt(:,:) = mbathy(:,:) - 1  
     326 
     327      ! Initialize / reset 
     328      ll_found(:,:) = .false. 
     329 
     330      IF ( rn_iso_frac - zepsilon > 0. ) THEN 
     331         ! Search for a uniform density/temperature region where adjacent levels           
     332         ! differ by less than rn_iso_frac * deltaT.                                       
     333         ! ik_iso is the index of the last level in the uniform layer   
     334         ! ll_found indicates whether the mixed layer depth can be found by interpolation  
     335         ik_iso(:,:)   = ik_ref(:,:)  
     336         DO jj = 1, nlcj  
     337            DO ji = 1, nlci  
     338!CDIR NOVECTOR  
     339               DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1  
     340                  IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN  
     341                     ik_iso(ji,jj)   = jk  
     342                     ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )  
     343                     EXIT  
     344                  END IF  
     345               END DO  
     346            END DO  
     347         END DO  
     348 
     349         ! Use linear interpolation to find depth of mixed layer base where possible  
     350         hmld_zint(:,:) = rn_zref  
     351         DO jj = 1, jpj  
     352            DO ji = 1, jpi  
     353               IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN  
     354                  zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )  
     355                  hmld_zint(ji,jj) = fsdept(ji,jj,ik_iso(ji,jj)) + zdz  
     356               END IF  
     357            END DO  
     358         END DO  
     359      END IF 
     360 
     361      ! If ll_found = .false. then calculate MLD using difference of zdelta_T     
     362      ! from the reference density/temperature  
     363  
     364! Prevent this section from working on land points  
     365      WHERE ( tmask(:,:,1) /= 1.0 )  
     366         ll_found = .true.  
     367      END WHERE  
     368  
     369      DO jk=1, jpk  
     370         ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:)   
     371      END DO  
     372  
     373! Set default value where interpolation cannot be used (ll_found=false)   
     374      DO jj = 1, jpj  
     375         DO ji = 1, jpi  
     376            IF ( .not. ll_found(ji,jj) )  hmld_zint(ji,jj) = fsdept(ji,jj,ikmt(ji,jj))  
     377         END DO  
     378      END DO  
     379 
     380      DO jj = 1, jpj  
     381         DO ji = 1, jpi  
     382!CDIR NOVECTOR  
     383            DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj)  
     384               IF ( ll_found(ji,jj) ) EXIT  
     385               IF ( ll_belowml(ji,jj,jk) ) THEN                 
     386                  zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) )  
     387                  zdT  = zT_b - zT(ji,jj,jk-1)                                       
     388                  zdz  = zdT / zdTdz(ji,jj,jk-1)                                        
     389                  hmld_zint(ji,jj) = fsdept(ji,jj,jk-1) + zdz  
     390                  EXIT                                                    
     391               END IF  
     392            END DO  
     393         END DO  
     394      END DO  
     395 
     396      hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1)  
     397      !   
     398      CALL wrk_dealloc( jpi, jpj, ikmt, ik_ref, ik_iso)  
     399      CALL wrk_dealloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 )  
     400      CALL wrk_dealloc( jpi,jpj, jpk, zT, zdTdz, zmoddT )  
     401      !  
     402   END SUBROUTINE zdf_mxl_zint_mld 
     403 
     404   SUBROUTINE zdf_mxl_zint_htc( kt ) 
     405      !!---------------------------------------------------------------------- 
     406      !!                  ***  ROUTINE zdf_mxl_zint_htc  *** 
     407      !!  
     408      !! ** Purpose :    
     409      !! 
     410      !! ** Method  :    
     411      !!---------------------------------------------------------------------- 
     412 
     413      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     414 
     415      INTEGER :: ji, jj, jk 
     416      INTEGER :: ikmax 
     417      REAL(wp) :: zc, zcoef 
     418      ! 
     419      INTEGER,  ALLOCATABLE, DIMENSION(:,:) ::   ilevel 
     420      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zthick_0, zthick 
     421 
     422      !!---------------------------------------------------------------------- 
     423 
     424      IF( .NOT. ALLOCATED(ilevel) ) THEN 
     425         ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), & 
     426         &         zthick(jpi,jpj), STAT=ji ) 
     427         IF( lk_mpp  )   CALL mpp_sum(ji) 
     428         IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' ) 
     429      ENDIF 
     430 
     431      ! Find last whole model T level above the MLD 
     432      ilevel(:,:)   = 0 
     433      zthick_0(:,:) = 0._wp 
     434 
     435      DO jk = 1, jpkm1   
     436         DO jj = 1, jpj 
     437            DO ji = 1, jpi                     
     438               zthick_0(ji,jj) = zthick_0(ji,jj) + fse3t(ji,jj,jk) 
     439               IF( zthick_0(ji,jj) < hmld_zint(ji,jj) )   ilevel(ji,jj) = jk 
     440            END DO 
     441         END DO 
     442         WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2) 
     443         WRITE(numout,*) 'fsdepw(jk+1 =',jk+1,') =',fsdepw(2,2,jk+1) 
     444      END DO 
     445 
     446      ! Surface boundary condition 
     447      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0._wp       ;   htc_mld(:,:) = 0._wp                                    
     448      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)    
     449      ENDIF 
     450 
     451      ! Deepest whole T level above the MLD 
     452      ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 ) 
     453 
     454      ! Integration down to last whole model T level 
     455      DO jk = 1, ikmax 
     456         DO jj = 1, jpj 
     457            DO ji = 1, jpi 
     458               zc = fse3t(ji,jj,jk) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1  )  )    ! 0 below ilevel 
     459               zthick(ji,jj) = zthick(ji,jj) + zc 
     460               htc_mld(ji,jj) = htc_mld(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 
     461            END DO 
     462         END DO 
     463      END DO 
     464 
     465      ! Subsequent partial T level 
     466      zthick(:,:) = hmld_zint(:,:) - zthick(:,:)   !   remaining thickness to reach MLD 
     467 
     468      DO jj = 1, jpj 
     469         DO ji = 1, jpi 
     470            htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem)  &  
     471      &                      * MIN( fse3t(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1) 
     472         END DO 
     473      END DO 
     474 
     475      WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2) 
     476 
     477      ! Convert to heat content 
     478      zcoef = rau0 * rcp 
     479      htc_mld(:,:) = zcoef * htc_mld(:,:) 
     480 
     481   END SUBROUTINE zdf_mxl_zint_htc 
     482 
     483   SUBROUTINE zdf_mxl_zint( kt ) 
     484      !!---------------------------------------------------------------------- 
     485      !!                  ***  ROUTINE zdf_mxl_zint  *** 
     486      !!  
     487      !! ** Purpose :    
     488      !! 
     489      !! ** Method  :    
     490      !!---------------------------------------------------------------------- 
     491 
     492      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     493 
     494      INTEGER :: ios 
     495      INTEGER :: jn 
     496 
     497      INTEGER :: nn_mld_diag = 0    ! number of diagnostics 
     498 
     499      INTEGER :: i_steps          ! no of timesteps per hour 
     500      INTEGER :: ierror             ! logical error message    
     501 
     502      REAL(wp) :: zdt             ! timestep variable   
     503 
     504      CHARACTER(len=1) :: cmld 
     505 
     506      TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 
     507      TYPE(MXL_ZINT), SAVE, DIMENSION(5) ::   mld_diags 
     508 
     509      NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 
     510 
     511      !!---------------------------------------------------------------------- 
     512       
     513      IF( kt == nit000 ) THEN 
     514         REWIND( numnam_ref )              ! Namelist namzdf_mldzint in reference namelist  
     515         READ  ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901) 
     516901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist', lwp ) 
     517 
     518         REWIND( numnam_cfg )              ! Namelist namzdf_mldzint in configuration namelist  
     519         READ  ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 ) 
     520902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist', lwp ) 
     521         IF(lwm) WRITE ( numond, namzdf_mldzint ) 
     522 
     523         IF( nn_mld_diag > 5 )   CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' ) 
     524 
     525         mld_diags(1) = sn_mld1 
     526         mld_diags(2) = sn_mld2 
     527         mld_diags(3) = sn_mld3 
     528         mld_diags(4) = sn_mld4 
     529         mld_diags(5) = sn_mld5 
     530 
     531         IF( nn_mld_diag > 0 ) THEN 
     532            WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================' 
     533            WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)' 
     534            DO jn = 1, nn_mld_diag 
     535               WRITE(numout,*) 'MLD criterion',jn,':' 
     536               WRITE(numout,*) '    nn_mld_type =', mld_diags(jn)%mld_type 
     537               WRITE(numout,*) '    rn_zref ='    , mld_diags(jn)%zref 
     538               WRITE(numout,*) '    rn_dT_crit =' , mld_diags(jn)%dT_crit 
     539               WRITE(numout,*) '    rn_iso_frac =', mld_diags(jn)%iso_frac 
     540            END DO 
     541            WRITE(numout,*) '====================================================================' 
     542         ENDIF 
     543      ENDIF 
     544 
     545      IF( nn_mld_diag > 0 ) THEN 
     546         DO jn = 1, nn_mld_diag 
     547            WRITE(cmld,'(I1)') jn 
     548            IF( iom_use( "mldzint_"//cmld ) .OR. iom_use( "mldhtc_"//cmld ) ) THEN 
     549               CALL zdf_mxl_zint_mld( mld_diags(jn) ) 
     550 
     551               IF( iom_use( "mldzint_"//cmld ) ) THEN 
     552                  CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) ) 
     553               ENDIF 
     554 
     555               IF( iom_use( "mldhtc_"//cmld ) )  THEN 
     556                  CALL zdf_mxl_zint_htc( kt ) 
     557                  CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:)   ) 
     558               ENDIF 
     559 
     560               IF( iom_use( "mldzint25h_"//cmld ) ) THEN 
     561                  IF( .NOT. mld_25h_write ) mld_25h_write = .TRUE. 
     562                  zdt = rdt 
     563                  IF( nacc == 1 ) zdt = rdtmin 
     564                  IF( MOD( 3600,INT(zdt) ) == 0 ) THEN 
     565                     i_steps = 3600/INT(zdt) 
     566                  ELSE 
     567                     CALL ctl_stop('STOP', 'zdf_mxl_zint 25h: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible') 
     568                  ENDIF 
     569                  IF( ( mld_25h_init ) .OR. ( kt == nit000 ) ) THEN 
     570                     i_cnt_25h = 1 
     571                     IF( .NOT. ALLOCATED(hmld_zint_25h) ) THEN 
     572                        ALLOCATE( hmld_zint_25h(jpi,jpj,nn_mld_diag), STAT=ierror ) 
     573                        IF( ierror > 0 )  CALL ctl_stop( 'zdf_mxl_zint 25h: unable to allocate hmld_zint_25h' )    
     574                     ENDIF 
     575                     hmld_zint_25h(:,:,jn) = hmld_zint(:,:) 
     576                  ENDIF 
     577                  IF( MOD( kt, i_steps ) == 0 .AND.  kt .NE. nn_it000 ) THEN 
     578                     hmld_zint_25h(:,:,jn) = hmld_zint_25h(:,:,jn) + hmld_zint(:,:) 
     579                  ENDIF 
     580                  IF( i_cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN 
     581                     CALL iom_put( "mldzint25h_"//cmld , hmld_zint_25h(:,:,jn) / 25._wp   ) 
     582                  ENDIF 
     583               ENDIF 
     584 
     585            ENDIF 
     586         END DO 
     587                   
     588         IF(  mld_25h_write  ) THEN 
     589            IF( ( MOD( kt, i_steps ) == 0 ) .OR.  mld_25h_init ) THEN 
     590               IF (lwp) THEN 
     591                  WRITE(numout,*) 'zdf_mxl_zint (25h) : Summed the following number of hourly values so far',i_cnt_25h 
     592          ENDIF 
     593               i_cnt_25h = i_cnt_25h + 1 
     594               IF( mld_25h_init ) mld_25h_init = .FALSE. 
     595            ENDIF 
     596            IF( i_cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN 
     597               i_cnt_25h = 1  
     598               DO jn = 1, nn_mld_diag  
     599                     hmld_zint_25h(:,:,jn) = hmld_zint(:,:)  
     600               ENDDO  
     601            ENDIF 
     602         ENDIF 
     603                   
     604      ENDIF 
     605 
     606   END SUBROUTINE zdf_mxl_zint 
    147607 
    148608   !!====================================================================== 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r10248 r10249  
    8585   USE stopar 
    8686   USE stopts 
     87   USE diatmb          ! Top,middle,bottom output 
     88   USE dia25h          ! 25h mean output 
     89   USE diaopfoam       ! FOAM operational output 
     90   USE diurnal_bulk    ! diurnal bulk SST  
    8791 
    8892   IMPLICIT NONE 
    … …  
    403407                            CALL  istate_init   ! ocean initial state (Dynamics and tracers) 
    404408 
     409      CALL diurnal_sst_bulk_init                ! diurnal sst 
     410      IF ( ln_diurnal ) CALL diurnal_sst_coolskin_init   ! cool skin   
     411 
    405412      IF( lk_tide       )   CALL    tide_init( nit000 )    ! Initialisation of the tidal harmonics 
    406413 
    … …  
    473480 
    474481      !                                     ! Assimilation increments 
    475       IF( lk_asminc     )   CALL asm_inc_init   ! Initialize assimilation increments 
     482      IF( lk_asminc ) THEN  
     483#if defined key_shelf  
     484         CALL  zdf_mxl_tref()     ! Initialization of hmld_tref 
     485#endif  
     486         CALL asm_inc_init     ! Initialize assimilation increments  
     487      ENDIF 
     488 
    476489      IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler 
     490                            CALL dia_tmb_init  ! TMB outputs 
     491                            CALL dia_25h_init  ! 25h mean  outputs 
     492                            CALL dia_diaopfoam_init  ! FOAM operational output 
    477493      ! 
    478494   END SUBROUTINE nemo_init 
    … …  
    630646      USE ldftra_oce, ONLY: ldftra_oce_alloc 
    631647      USE trc_oce   , ONLY: trc_oce_alloc 
     648      USE diainsitutem, ONLY: insitu_tem_alloc 
    632649#if defined key_diadct  
    633650      USE diadct    , ONLY: diadct_alloc  
    … …  
    646663      ierr = ierr + ldftra_oce_alloc()          ! ocean lateral  physics : tracers 
    647664      ierr = ierr + zdf_oce_alloc   ()          ! ocean vertical physics 
     665      ierr = ierr + insitu_tem_alloc() 
    648666      ! 
    649667      ierr = ierr + trc_oce_alloc   ()          ! shared TRC / TRA arrays 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/step.F90

    r10248 r10249  
    225225      ENDIF 
    226226 
     227      ! Cool skin 
     228      IF ( ln_diurnal ) THEN   
     229         IF ( ln_blk_core ) THEN 
     230            CALL diurnal_sst_coolskin_step( &   
     231                    qns(:,:)+(rn_abs*qsr(:,:)), taum, rhop(:,:,1), rdt)  
     232         ELSE 
     233            CALL diurnal_sst_coolskin_step( &   
     234                    qns, taum, rhop(:,:,1), rdt)  
     235         ENDIF 
     236      ENDIF 
     237 
    227238      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    228239      ! diagnostics and outputs             (ua, va, tsa used as workspace) 
    … …  
    238249      IF( ln_crs     )      CALL crs_fld( kstp )         ! ocean model: online field coarsening & output 
    239250 
     251      !Diurnal warm layer model         
     252      IF ( ln_diurnal ) THEN 
     253         IF ( ln_blk_core ) THEN 
     254            IF( kstp == nit000 )THEN   
     255               CALL diurnal_sst_takaya_step( &   
     256               &    qsr(:,:)-(rn_abs*qsr(:,:)), qns(:,:)+(rn_abs*qsr(:,:)), & 
     257               &    taum, rhop(:,:,1), & 
     258               &    rdt, ld_calcfrac = .TRUE.)   
     259            ELSE   
     260               CALL diurnal_sst_takaya_step( &   
     261               &    qsr(:,:)-(rn_abs*qsr(:,:)), qns(:,:)+(rn_abs*qsr(:,:)), & 
     262               &    taum, rhop(:,:,1), rdt )   
     263            ENDIF  
     264         ELSE 
     265            IF( kstp == nit000 )THEN   
     266               CALL diurnal_sst_takaya_step( &   
     267               &    qsr, qns, taum, rhop(:,:,1), & 
     268               &    rdt, ld_calcfrac = .TRUE.)   
     269            ELSE   
     270               CALL diurnal_sst_takaya_step( &   
     271               &    qsr, qns, taum, rhop(:,:,1), rdt )   
     272            ENDIF  
     273         ENDIF 
     274      ENDIF 
     275 
    240276#if defined key_top 
    241277      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    … …  
    255291                             CALL tra_sbc    ( kstp )       ! surface boundary condition 
    256292      IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr 
     293      IF( ln_tradwl      )   CALL tra_dwl    ( kstp )       ! Polcoms Style Short Wave Radiation  
    257294      IF( ln_trabbc      )   CALL tra_bbc    ( kstp )       ! bottom heat flux 
    258295      IF( lk_trabbl      )   CALL tra_bbl    ( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme 
    … …  
    299336      ! Dynamics                                    (tsa used as workspace) 
    300337      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     338 
     339      IF( ln_bkgwri )        CALL asm_bkg_wri( kstp )     ! output background fields 
     340 
    301341      IF( lk_dynspg_ts   )  THEN 
    302342                                                             ! revert to previously computed momentum tendencies 
    … …  
    317357        IF(  lk_asminc .AND. ln_asmiau .AND. & 
    318358           & ln_dyninc      )  CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment 
    319         IF( ln_bkgwri )        CALL asm_bkg_wri( kstp )     ! output background fields 
    320359        IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified) 
    321360        IF( lk_bdy          )  CALL bdy_dyn3d_dmp(kstp )    ! bdy damping trends 
    … …  
    336375                               CALL ssh_swp( kstp )         ! swap of sea surface height 
    337376      IF( lk_vvl           )   CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors 
    338       ! 
    339       IF( lrst_oce         )   CALL rst_write( kstp )       ! write output ocean restart file 
    340377 
    341378#if defined key_agrif 
    … …  
    361398                               CALL dia_wri_state( 'output.abort', kstp ) 
    362399      ENDIF 
     400      IF( ln_harm_ana_store   )   CALL harm_ana( kstp )        ! Harmonic analysis of tides  
    363401      IF( kstp == nit000   )   THEN 
    364402                 CALL iom_close( numror )     ! close input  ocean restart file 
    … …  
    366404         IF( lwm.AND.numoni /= -1 ) CALL FLUSH    ( numoni )     ! flush output namelist ice 
    367405      ENDIF 
     406      IF( lrst_oce         )   CALL rst_write    ( kstp )   ! write output ocean restart file 
    368407 
    369408      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r10248 r10249  
    2525   USE sbcrnf           ! surface boundary condition: runoff variables 
    2626   USE sbccpl           ! surface boundary condition: coupled formulation (call send at end of step) 
     27   USE sbcflx           ! surface boundary condition: Fluxes 
    2728   USE sbc_oce          ! surface boundary condition: ocean 
    2829   USE sbctide          ! Tide initialisation 
    … …  
    3031 
    3132   USE traqsr           ! solar radiation penetration      (tra_qsr routine) 
     33   USE tradwl           ! POLCOMS style solar radiation    (tra_dwl routine)  
    3234   USE trasbc           ! surface boundary condition       (tra_sbc routine) 
    3335   USE trabbc           ! bottom boundary condition        (tra_bbc routine) 
    … …  
    106108   USE prtctl           ! Print control                    (prt_ctl routine) 
    107109 
     110   USE harmonic_analysis ! harmonic analysis of tides (harm_ana routine)  
     111   USE bdytides          ! harmonic analysis of tides (harm_ana routine)  
    108112   USE diaobs           ! Observation operator 
     113 
     114   USE diurnal_bulk    ! diurnal SST bulk routines  (diurnal_sst_takaya routine)   
     115   USE cool_skin       ! diurnal cool skin correction (diurnal_sst_coolskin routine) 
    109116 
    110117   USE timing           ! Timing 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90

    r10246 r10249  
    2727   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   etot3         !: light absortion coefficient 
    2828   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   facvol        !: volume for degraded regions 
     29   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   rlambda2      !: Lambda2 for downwell version of Short wave Radiation 
     30   REAL(wp), PUBLIC                                      ::   rlambda       !: Lambda  for downwell version of Short wave Radiation 
    2931 
    3032#if defined key_top  
    … …  
    7880      !!                  ***  trc_oce_alloc  *** 
    7981      !!---------------------------------------------------------------------- 
    80       INTEGER ::   ierr(2)        ! Local variables 
     82      INTEGER ::   ierr(3)        ! Local variables 
    8183      !!---------------------------------------------------------------------- 
    8284      ierr(:) = 0 
    8385                     ALLOCATE( etot3 (jpi,jpj,jpk), STAT=ierr(1) ) 
    8486      IF( lk_degrad) ALLOCATE( facvol(jpi,jpj,jpk), STAT=ierr(2) ) 
     87                    ALLOCATE( rlambda2(jpi,jpj),   STAT=ierr(3) ) 
    8588      trc_oce_alloc  = MAXVAL( ierr ) 
    8689      ! 
    8790      IF( trc_oce_alloc /= 0 )   CALL ctl_warn('trc_oce_alloc: failed to allocate etot3 array') 
     91      IF( trc_oce_alloc /= 0 )   CALL ctl_warn('trc_oce_alloc: failed to allocate etot3, facvol or rlambda2 array') 
    8892   END FUNCTION trc_oce_alloc 
    8993 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcsms_my_trc.F90

    r10246 r10249  
    1818   USE trd_oce 
    1919   USE trdtrc 
     20   USE trcbc, only : trc_bc_read 
    2021 
    2122   IMPLICIT NONE 
    … …  
    5556 
    5657      IF( l_trdtrc )  CALL wrk_alloc( jpi, jpj, jpk, ztrmyt ) 
     58 
     59      CALL trc_bc_read  ( kt )       ! tracers: surface and lateral Boundary Conditions 
    5760 
    5861      IF( l_trdtrc ) THEN      ! Save the trends in the ixed layer 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/TOP_SRC/MY_TRC/trcwri_my_trc.F90

    r10246 r10249  
    1919 
    2020   PUBLIC trc_wri_my_trc  
     21#if defined key_tracer_budget 
     22   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: trb_temp ! slwa 
     23#endif 
     24 
    2125 
    2226#  include "top_substitute.h90" 
    2327CONTAINS 
    2428 
     29#if defined key_tracer_budget 
     30   SUBROUTINE trc_wri_my_trc (kt, fl) ! slwa 
     31#else 
    2532   SUBROUTINE trc_wri_my_trc 
     33#endif 
    2634      !!--------------------------------------------------------------------- 
    2735      !!                     ***  ROUTINE trc_wri_trc  *** 
    … …  
    2937      !! ** Purpose :   output passive tracers fields  
    3038      !!--------------------------------------------------------------------- 
     39#if defined key_tracer_budget 
     40      INTEGER, INTENT( in ), OPTIONAL     :: fl  
     41      INTEGER, INTENT( in )               :: kt 
     42      REAL(wp), DIMENSION(jpi,jpj,jpk)    :: trpool !tracer pool temporary output 
     43#else 
     44      INTEGER, INTENT( in )               :: kt 
     45#endif 
    3146      CHARACTER (len=20)   :: cltra 
    32       INTEGER              :: jn 
     47      INTEGER              :: jn,jk ! JC TODO jk defined here but may not be used 
    3348      !!--------------------------------------------------------------------- 
    3449  
    3550      ! write the tracer concentrations in the file 
    3651      ! --------------------------------------- 
     52 
     53 
     54#if defined key_tracer_budget 
     55      IF( PRESENT(fl)) THEN 
     56! depth integrated 
     57! for strict budgetting write this out at end of timestep as an average between 'now' and 'after' at kt 
     58         DO jn = jp_myt0, jp_myt1  
     59          IF(ln_trdtrc (jn))THEN 
     60            trpool(:,:,:) = 0.5 * ( trn(:,:,:,jn) * fse3t_a(:,:,:) +  & 
     61                                        trb_temp(:,:,:,jn) * fse3t(:,:,:) ) 
     62            cltra = TRIM( ctrcnm(jn) )//"e3t"     ! depth integrated output 
     63            IF( kt == nittrc000 ) write(6,*)'output pool ',cltra 
     64            DO jk = 1, jpk 
     65               trpool(:,:,jk) = trpool(:,:,jk) 
     66            END DO 
     67            CALL iom_put( cltra, trpool) 
     68 
     69          END IF 
     70         END DO 
     71 
     72      ELSE 
     73 
     74         IF( kt == nittrc000 ) THEN 
     75           ALLOCATE(trb_temp(jpi,jpj,jpk,jp_my_trc))  ! slwa 
     76         ENDIF 
     77         trb_temp(:,:,:,:)=trn(:,:,:,:) ! slwa save for tracer budget (unfiltered trn) 
     78 
     79 
     80      END IF 
     81#else 
    3782      DO jn = jp_myt0, jp_myt1 
    3883         cltra = TRIM( ctrcnm(jn) )                  ! short title for tracer 
    3984         CALL iom_put( cltra, trn(:,:,:,jn) ) 
    4085      END DO 
     86#endif 
    4187      ! 
    4288   END SUBROUTINE trc_wri_my_trc 
    … …  
    4894   PUBLIC trc_wri_my_trc 
    4995CONTAINS 
    50    SUBROUTINE trc_wri_my_trc                     ! Empty routine   
     96#if defined key_tracer_budget 
     97   SUBROUTINE trc_wri_my_trc (kt, fl) ! slwa 
     98      INTEGER, INTENT( in ), OPTIONAL     :: fl  
     99      INTEGER, INTENT( in )               :: kt 
     100#else 
     101   ! JC TODO Subroutine arguments (kt) inconsistent with earlier definition 
     102   SUBROUTINE trc_wri_my_trc (kt) 
     103      INTEGER, INTENT( in )               :: kt 
     104#endif 
    51105   END SUBROUTINE trc_wri_my_trc 
    52106#endif 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90

    r10246 r10249  
    5656      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    5757      !! 
    58       INTEGER            :: jn 
     58      INTEGER            :: jn, jk 
    5959      CHARACTER (len=22) :: charout 
    6060      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrtrd 
    … …  
    105105        DO jn = 1, jptra 
    106106           ztrtrd(:,:,:,jn) = tra(:,:,:,jn) - ztrtrd(:,:,:,jn) 
     107#if defined key_tracer_budget 
     108           DO jk = 1, jpkm1 
     109             ztrtrd(:,:,jk,jn) = ztrtrd(:,:,jk,jn) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk)  ! slwa 
     110           END DO 
     111#endif 
    107112           CALL trd_tra( kt, 'TRC', jn, jptra_ldf, ztrtrd(:,:,:,jn) ) 
    108113        END DO 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90

    r10248 r10249  
    3333   USE trdtra 
    3434   USE tranxt 
     35   USE trcbdy          ! BDY open boundaries 
     36   USE bdy_par, only: lk_bdy 
     37   USE iom 
    3538# if defined key_agrif 
    3639   USE agrif_top_interp 
    … …  
    9396      CHARACTER (len=22) :: charout 
    9497      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::  ztrdt  
     98#if defined key_tracer_budget 
     99      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  ztrdt_m1 ! slwa 
     100#endif 
    95101      !!---------------------------------------------------------------------- 
    96102      ! 
    … …  
    101107         WRITE(numout,*) 'trc_nxt : time stepping on passive tracers' 
    102108      ENDIF 
     109#if defined key_tracer_budget 
     110      IF( kt == nittrc000 .AND. l_trdtrc ) THEN 
     111         ALLOCATE( ztrdt_m1(jpi,jpj,jpk,jptra) )  ! slwa 
     112         IF( ln_rsttr .AND.    &                     ! Restart: read in restart  file 
     113            iom_varid( numrtr, 'atf_trend_'//TRIM(ctrcnm(1)), ldstop = .FALSE. ) > 0 ) THEN 
     114            IF(lwp) WRITE(numout,*) '          nittrc000-nn_dttrc ATF tracer trend read in the restart file' 
     115            DO jn = 1, jptra 
     116               CALL iom_get( numrtr, jpdom_autoglo, 'atf_trend_'//TRIM(ctrcnm(jn)), ztrdt_m1(:,:,:,jn) )   ! before tracer trend for atf 
     117            END DO 
     118         ELSE           
     119           ztrdt_m1=0.0 
     120         ENDIF 
     121      ENDIF 
     122#endif 
     123 
    103124 
    104125#if defined key_agrif 
    … …  
    111132 
    112133 
    113 #if defined key_bdy 
    114 !!      CALL bdy_trc( kt )               ! BDY open boundaries 
    115 #endif 
     134      IF( lk_bdy )  CALL trc_bdy( kt )               ! BDY open boundaries 
    116135 
    117136 
    … …  
    149168               zfact = 1.e0 / r2dt(jk)   
    150169               ztrdt(:,:,jk,jn) = ( trb(:,:,jk,jn) - ztrdt(:,:,jk,jn) ) * zfact  
    151                CALL trd_tra( kt, 'TRC', jn, jptra_atf, ztrdt ) 
     170!slwa          CALL trd_tra( kt, 'TRC', jn, jptra_atf, ztrdt ) 
     171#if defined key_tracer_budget 
     172               ztrdt(:,:,jk,jn) = ztrdt(:,:,jk,jn) * e1t(:,:) * e2t(:,:) * e3t_n(:,:,jk)  ! slwa vvl 
     173               !ztrdt(:,:,jk,jn) = ztrdt(:,:,jk,jn) * e1t(:,:) * e2t(:,:) * e3t_0(:,:,jk)  ! slwa CHANGE for vvl 
     174#endif 
    152175            END DO 
     176#if defined key_tracer_budget 
     177! slwa budget code 
     178              CALL trd_tra( kt, 'TRC', jn, jptra_atf, ztrdt_m1(:,:,:,jn) ) 
     179#else 
     180              CALL trd_tra( kt, 'TRC', jn, jptra_atf, ztrdt(:,:,:,jn) ) 
     181#endif 
    153182         END DO 
     183#if defined key_tracer_budget 
     184        ztrdt_m1(:,:,:,:) = ztrdt(:,:,:,:)    ! need previous time step for budget slwa 
     185#endif 
    154186         CALL wrk_dealloc( jpi, jpj, jpk, jptra, ztrdt )  
    155187      END IF 
     188 
     189#if defined key_tracer_budget 
     190      !                                           Write in the tracer restart file 
     191      !                                          ******************************* 
     192      IF( lrst_trc ) THEN 
     193         IF(lwp) WRITE(numout,*) 
     194         IF(lwp) WRITE(numout,*) 'trc : ATF trend at last time step for tracer budget written in tracer restart file ',   & 
     195            &                    'at it= ', kt,' date= ', ndastp 
     196         IF(lwp) WRITE(numout,*) '~~~~' 
     197         DO jn = 1, jptra 
     198            CALL iom_rstput( kt, nitrst, numrtw, 'atf_trend_'//TRIM(ctrcnm(jn)), ztrdt_m1(:,:,:,jn) ) 
     199         END DO 
     200      ENDIF 
     201#endif 
     202 
    156203      ! 
    157204      IF(ln_ctl)   THEN  ! print mean trends (used for debugging) 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/TOP_SRC/TRP/trcrad.F90

    r10246 r10249  
    1818   USE trdtra 
    1919   USE prtctl_trc          ! Print control for debbuging 
     20#if defined key_tracer_budget 
     21   USE iom 
     22#endif 
    2023 
    2124   IMPLICIT NONE 
    … …  
    5154      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index       
    5255      CHARACTER (len=22) :: charout 
     56      ! +++>>> FABM 
     57      INTEGER :: jn 
     58      ! FABM <<<+++ 
    5359      !!---------------------------------------------------------------------- 
    5460      ! 
    … …  
    6571      IF( lk_pisces  )   CALL trc_rad_sms( kt, trb, trn, jp_pcs0 , jp_pcs1, cpreserv='Y' )  ! PISCES model 
    6672      IF( lk_my_trc  )   CALL trc_rad_sms( kt, trb, trn, jp_myt0 , jp_myt1               )  ! MY_TRC model 
    67  
     73      ! +++>>> FABM 
     74      IF( lk_fabm  )   THEN 
     75        DO jn=1,jp_fabm ! state variable loop 
     76          IF (lk_rad_fabm(jn)) THEN 
     77           CALL trc_rad_sms( kt, trb, trn, jn+jp_fabm_m1 , jn+jp_fabm_m1 ) 
     78          ENDIF 
     79        END DO 
     80      END IF 
     81      ! FABM <<<+++ 
    6882      ! 
    6983      IF(ln_ctl) THEN      ! print mean trends (used for debugging) 
    … …  
    110124      REAL(wp) :: zcoef, ztrcorn, ztrmasn   !    "         " 
    111125      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrtrdb, ztrtrdn   ! workspace arrays 
     126#if defined key_tracer_budget 
     127      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  ztrtrdb_m1 ! slwa  
     128#endif 
    112129      REAL(wp) :: zs2rdt 
    113130      LOGICAL ::   lldebug = .FALSE. 
    … …  
    116133  
    117134      IF( l_trdtrc )  CALL wrk_alloc( jpi, jpj, jpk, ztrtrdb, ztrtrdn ) 
     135#if defined key_tracer_budget 
     136      IF( kt == nittrc000 .AND. l_trdtrc) THEN 
     137         IF (.not. ALLOCATED(ztrtrdb_m1)) ALLOCATE( ztrtrdb_m1(jpi,jpj,jpk,jptra) )  ! slwa 
     138            DO jn = jp_sms0, jp_sms1 
     139               IF( ln_rsttr .AND.    &                     ! Restart: read in restart  file 
     140                  iom_varid( numrtr, 'rdb_trend_'//TRIM(ctrcnm(jn)), ldstop = .FALSE. ) > 0 ) THEN 
     141                  IF(lwp) WRITE(numout,*) '          nittrc000-nn_dttrc RDB tracer trend read for',TRIM(ctrcnm(jn)) 
     142                  CALL iom_get( numrtr, jpdom_autoglo, 'rdb_trend_'//TRIM(ctrcnm(jn)), ztrtrdb_m1(:,:,:,jn) )   ! before tracer trend for rdb 
     143               ELSE 
     144                  IF(lwp) WRITE(numout,*) '          no nittrc000-nn_dttrc RDB tracer trend for',TRIM(ctrcnm(jn)),', setting to 0.' 
     145                  ztrtrdb_m1(:,:,:,jn)=0.0 
     146               ENDIF 
     147            END DO 
     148      ENDIF 
     149#endif 
    118150       
    119151      IF( PRESENT( cpreserv )  ) THEN   !  total tracer concentration is preserved  
    … …  
    156188               ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt 
    157189               ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt  
     190#if defined key_tracer_budget 
     191! slwa budget code 
     192               DO jk = 1, jpkm1 
     193                  ztrtrdb(:,:,jk) = ztrtrdb(:,:,jk) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk) 
     194                  ztrtrdn(:,:,jk) = ztrtrdn(:,:,jk) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk) 
     195               END DO 
     196               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb_m1(:,:,:,jn) ) 
     197               ztrtrdb_m1(:,:,:,jn)=ztrtrdb(:,:,:) 
     198#else 
    158199               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb )       ! Asselin-like trend handling 
     200#endif 
    159201               CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn )       ! standard     trend handling 
    160202              ! 
    … …  
    187229               ztrtrdb(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrdb(:,:,:) ) * zs2rdt 
    188230               ztrtrdn(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrdn(:,:,:) ) * zs2rdt  
     231#if defined key_tracer_budget 
     232! slwa budget code 
     233               DO jk = 1, jpkm1 
     234                  ztrtrdb(:,:,jk) = ztrtrdb(:,:,jk) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk) 
     235                  ztrtrdn(:,:,jk) = ztrtrdn(:,:,jk) * e1t(:,:) * e2t(:,:) * fse3t(:,:,jk) 
     236               END DO 
     237               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb_m1(:,:,:,jn) ) 
     238               ztrtrdb_m1(:,:,:,jn)=ztrtrdb(:,:,:) 
     239#else 
    189240               CALL trd_tra( kt, 'TRC', jn, jptra_radb, ztrtrdb )       ! Asselin-like trend handling 
     241#endif 
    190242               CALL trd_tra( kt, 'TRC', jn, jptra_radn, ztrtrdn )       ! standard     trend handling 
    191243              ! 
    … …  
    195247 
    196248      ENDIF 
     249 
     250#if defined key_tracer_budget 
     251      !                                           Write in the tracer restart file 
     252      !                                          ******************************* 
     253      IF( lrst_trc ) THEN 
     254         IF(lwp) WRITE(numout,*) 
     255         IF(lwp) WRITE(numout,*) 'trc : RDB trend at last time step for tracer budget written in tracer restart file ',   & 
     256            &                    'at it= ', kt,' date= ', ndastp 
     257         IF(lwp) WRITE(numout,*) '~~~~' 
     258         DO jn = jp_sms0, jp_sms1 
     259            CALL iom_rstput( kt, nitrst, numrtw, 'rdb_trend_'//TRIM(ctrcnm(jn)), ztrtrdb_m1(:,:,:,jn) ) 
     260         END DO 
     261      ENDIF 
     262#endif 
    197263 
    198264      IF( l_trdtrc )  CALL wrk_dealloc( jpi, jpj, jpk, ztrtrdb, ztrtrdn ) 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90

    r10246 r10249  
    113113           sbc_trc_b(:,:,:) = 0._wp 
    114114         ENDIF 
     115         sbc_trc(:,:,:) = 0._wp    !slwa initialise for vvl 
    115116      ELSE                                         ! Swap of forcing fields 
    116117         IF( ln_top_euler ) THEN 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90

    r10246 r10249  
    2727   USE trcsbc          ! surface boundary condition          (trc_sbc routine) 
    2828   USE zpshde          ! partial step: hor. derivative       (zps_hde routine) 
     29   USE trcbdy          ! BDY open boundaries 
     30   USE bdy_par, only: lk_bdy 
    2931 
    3032#if defined key_agrif 
    … …  
    6870         IF( ln_trcdmp )        CALL trc_dmp( kstp )            ! internal damping trends 
    6971         IF( ln_trcdmp_clo )    CALL trc_dmp_clo( kstp )        ! internal damping trends on closed seas only 
     72         IF( lk_bdy )           CALL trc_bdy_dmp( kstp )        ! BDY damping trends 
    7073                                CALL trc_adv( kstp )            ! horizontal & vertical advection  
    7174                                CALL trc_ldf( kstp )            ! lateral mixing 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/TOP_SRC/TRP/trdtrc.F90

    r10246 r10249  
    102102         END SELECT 
    103103                                          cltra = TRIM(cltra)//TRIM(ctrcnm(kjn)) 
     104         ! +++>>>FABM 
     105#if defined key_tracer_budget 
     106! for outputting depth integrated 
     107         SELECT CASE( ktrd ) 
     108         CASE( jptra_xad, jptra_yad, jptra_zad  )  
     109           cltra = TRIM(cltra)//"_e3t" 
     110         END SELECT 
     111#endif 
     112         ! FABM <<<+++ 
    104113                                          CALL iom_put( cltra,  ptrtrd(:,:,:) ) 
    105114         ! 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/TOP_SRC/TRP/trdtrc_oce.F90

    r10246 r10249  
    2222   CHARACTER(len=50) ::  cn_trdrst_trc_in     !: suffix of pass. tracer restart name (input) 
    2323   CHARACTER(len=50) ::  cn_trdrst_trc_out    !: suffix of pass. tracer restart name (output) 
    24    LOGICAL, DIMENSION(jptra) ::   ln_trdtrc   !: large trends diagnostic to write or not (namelist) 
     24   ! --->>> FABM 
     25   ! LOGICAL, DIMENSION(jptra) ::   ln_trdtrc   !: large trends diagnostic to write or not (namelist) 
     26   ! FABM <<<--- 
     27   ! +++>>> FABM 
     28   LOGICAL, DIMENSION(jpmaxtrc) ::   ln_trdtrc   !: large trends diagnostic to write or not (namelist) 
     29   ! FABM <<<+++ 
    2530 
    2631# if defined key_trdtrc && defined key_iomput 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/TOP_SRC/par_trc.F90

    r10246 r10249  
    1515   USE par_cfc       ! CFC 11 and 12 tracers 
    1616   USE par_my_trc    ! user defined passive tracers 
     17   ! +++>>> FABM 
     18   USE par_fabm      ! FABM 
     19   ! FABM <<<+++ 
    1720 
    1821   IMPLICIT NONE 
    … …  
    2427   ! Passive tracers : Total size 
    2528   ! ---------------               ! total number of passive tracers, of 2d and 3d output and trend arrays 
    26    INTEGER, PUBLIC,  PARAMETER ::   jptra    =  jp_pisces     + jp_cfc     + jp_c14b    + jp_my_trc 
    27    INTEGER, PUBLIC,  PARAMETER ::   jpdia2d  =  jp_pisces_2d  + jp_cfc_2d  + jp_c14b_2d + jp_my_trc_2d 
    28    INTEGER, PUBLIC,  PARAMETER ::   jpdia3d  =  jp_pisces_3d  + jp_cfc_3d  + jp_c14b_3d + jp_my_trc_3d 
     29   ! --->>> FABM 
     30   ! INTEGER, PUBLIC,  PARAMETER ::   jptra    =  jp_pisces     + jp_cfc     + jp_c14b    + jp_my_trc 
     31   ! INTEGER, PUBLIC,  PARAMETER ::   jpdia2d  =  jp_pisces_2d  + jp_cfc_2d  + jp_c14b_2d + jp_my_trc_2d 
     32   ! INTEGER, PUBLIC,  PARAMETER ::   jpdia3d  =  jp_pisces_3d  + jp_cfc_3d  + jp_c14b_3d + jp_my_trc_3d 
     33   ! FABM <<<--- 
     34   ! +++>>> FABM 
     35   INTEGER, PUBLIC  ::   jptra    =  jp_pisces     + jp_cfc     + jp_c14b    + jp_my_trc 
     36   INTEGER, PUBLIC  ::   jpdia2d  =  jp_pisces_2d  + jp_cfc_2d  + jp_c14b_2d + jp_my_trc_2d 
     37   INTEGER, PUBLIC  ::   jpdia3d  =  jp_pisces_3d  + jp_cfc_3d  + jp_c14b_3d + jp_my_trc_3d 
     38   ! FABM <<<+++ 
    2939   !                     ! total number of sms diagnostic arrays 
    30    INTEGER, PUBLIC,  PARAMETER ::   jpdiabio =  jp_pisces_trd + jp_cfc_trd + jp_c14b_trd + jp_my_trc_trd 
     40   ! --->>> FABM 
     41   ! INTEGER, PUBLIC,  PARAMETER ::   jpdiabio =  jp_pisces_trd + jp_cfc_trd + jp_c14b_trd + jp_my_trc_trd 
     42   ! FABM <<<--- 
     43   ! +++>>> FABM 
     44   INTEGER, PUBLIC  ::   jpdiabio =  jp_pisces_trd + jp_cfc_trd + jp_c14b_trd + jp_my_trc_trd 
     45   ! FABM <<<+++ 
    3146    
    3247   !  1D configuration ("key_c1d") 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/TOP_SRC/trc.F90

    r10246 r10249  
    1414   USE par_oce 
    1515   USE par_trc 
     16#if defined key_bdy 
     17   USE bdy_oce, only: nb_bdy, OBC_DATA 
     18#endif 
    1619    
    1720   IMPLICIT NONE 
    … …  
    8083   END TYPE 
    8184 
    82    REAL(wp), DIMENSION(jptra), PUBLIC         :: trc_ice_ratio, & ! ice-ocean tracer ratio 
     85   ! --->>> FABM  
     86   !REAL(wp), DIMENSION(jptra), PUBLIC         :: trc_ice_ratio, & ! ice-ocean tracer ratio 
     87   !                                              trc_ice_prescr   ! prescribed ice trc cc 
     88   !CHARACTER(len=2), DIMENSION(jptra), PUBLIC :: cn_trc_o ! choice of ocean tracer cc 
     89   ! FABM <<<--- 
     90   ! +++>>> FABM  
     91   REAL(wp), DIMENSION(jpmaxtrc), PUBLIC         :: trc_ice_ratio, & ! ice-ocean tracer ratio 
    8392                                                 trc_ice_prescr   ! prescribed ice trc cc 
    84    CHARACTER(len=2), DIMENSION(jptra), PUBLIC :: cn_trc_o ! choice of ocean tracer cc 
     93   CHARACTER(len=2), DIMENSION(jpmaxtrc), PUBLIC :: cn_trc_o ! choice of ocean tracer cc 
     94   ! FABM <<<+++ 
    8595 
    8696   !! information for outputs 
    … …  
    90100       CHARACTER(len = 80)  :: cllname  !: long name 
    91101       CHARACTER(len = 20)  :: clunit   !: unit 
    92        LOGICAL              :: llinit   !: read in a file or not 
    93        LOGICAL              :: llsave   !: save the tracer or not 
     102! --->>> FABM 
     103!       LOGICAL              :: llinit   !: read in a file or not 
     104!!#if defined  key_my_trc 
     105!       LOGICAL              :: llsbc   !: read in a file or not 
     106!       LOGICAL              :: llcbc   !: read in a file or not 
     107!       LOGICAL              :: llobc   !: read in a file or not 
     108!#endif 
     109!       LOGICAL              :: llsave   !: save the tracer or not 
     110! FABM <<<--- 
     111! +++ FABM 
     112       LOGICAL              :: llinit=.FALSE.   !: read in a file or not 
     113#if defined  key_fabm 
     114       LOGICAL              :: llsbc=.FALSE.   !: read in a file or not 
     115       LOGICAL              :: llcbc=.FALSE.   !: read in a file or not 
     116       LOGICAL              :: llobc=.FALSE.   !: read in a file or not 
     117#endif 
     118       LOGICAL              :: llsave=.FALSE.   !: save the tracer or not 
     119! FABM <<<+++ 
    94120   END TYPE PTRACER 
    95121   CHARACTER(len = 20), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)    ::  ctrcnm         !: tracer name  
    … …  
    191217# endif 
    192218   ! 
     219#if defined key_bdy 
     220   CHARACTER(len=20), PUBLIC, ALLOCATABLE,  SAVE,  DIMENSION(:)   ::  cn_trc_dflt          ! Default OBC condition for all tracers 
     221   CHARACTER(len=20), PUBLIC, ALLOCATABLE,  SAVE,  DIMENSION(:)   ::  cn_trc               ! Choice of boundary condition for tracers 
     222   INTEGER,           PUBLIC, ALLOCATABLE,  SAVE,  DIMENSION(:)   ::  nn_trcdmp_bdy        !: =T Tracer damping 
     223   ! External data structure of BDY for TOP. Available elements: cn_obc, ll_trc, trcnow, dmp 
     224   TYPE(OBC_DATA),    PUBLIC, ALLOCATABLE, DIMENSION(:,:), TARGET ::  trcdta_bdy           !: bdy external data (local process) 
     225#endif 
    193226 
    194227   !!---------------------------------------------------------------------- 
    … …  
    213246         &      cvol(jpi,jpj,jpk)     , rdttrc(jpk)           , trai(jptra)           ,       & 
    214247         &      ctrcnm(jptra)         , ctrcln(jptra)         , ctrcun(jptra)         ,       &  
     248! --->>> FABM 
     249!!#if defined key_my_trc 
     250! FABM <<<--- 
     251! +++>>> FABM 
     252#if defined key_fabm 
     253! FABM <<<+++ 
     254         &      ln_trc_sbc(jptra)     , ln_trc_cbc(jptra)     , ln_trc_obc(jptra)     ,       & 
     255#endif 
     256#if defined key_bdy 
     257         &      cn_trc_dflt(nb_bdy)   , cn_trc(nb_bdy)        , nn_trcdmp_bdy(nb_bdy) ,       & 
     258         &      trcdta_bdy(jptra,nb_bdy)                                              ,       & 
     259#endif 
    215260         &      ln_trc_ini(jptra)     , ln_trc_wri(jptra)     , qsr_mean(jpi,jpj)     ,  STAT = trc_alloc  )   
    216261 

  • branches/UKMO/dev_r5518_AMM15_package/NEMOGCM/NEMO/TOP_SRC/trcbc.F90

    r10246 r10249  
    44   !! TOP :  module for passive tracer boundary conditions 
    55   !!===================================================================== 
    6    !!---------------------------------------------------------------------- 
    7 #if  defined key_top  
     6   !! History :  3.5 !  2014-04  (M. Vichi, T. Lovato)  Original 
     7   !!            3.6 !  2015-03  (T . Lovato) Revision and BDY support 
     8   !!---------------------------------------------------------------------- 
     9#if defined key_top 
    810   !!---------------------------------------------------------------------- 
    911   !!   'key_top'                                                TOP model  
    1012   !!---------------------------------------------------------------------- 
    11    !!   trc_dta    : read and time interpolated passive tracer data 
     13   !!   trc_bc       : read and time interpolated tracer Boundary Conditions 
    1214   !!---------------------------------------------------------------------- 
    1315   USE par_trc       !  passive tracers parameters 
    … …  
    1719   USE lib_mpp       !  MPP library 
    1820   USE fldread       !  read input fields 
     21#if defined key_bdy 
     22   USE bdy_oce, only: nb_bdy , idx_bdy, ln_coords_file, rn_time_dmp, rn_time_dmp_out 
     23#endif 
    1924 
    2025   IMPLICIT NONE 
    … …  
    2429   PUBLIC   trc_bc_read    ! called in trcstp.F90 or within 
    2530 
    26    INTEGER  , SAVE, PUBLIC                             :: nb_trcobc   ! number of tracers with open BC 
    27    INTEGER  , SAVE, PUBLIC                             :: nb_trcsbc   ! number of tracers with surface BC 
    28    INTEGER  , SAVE, PUBLIC                             :: nb_trccbc   ! number of tracers with coastal BC 
     31   INTEGER  , SAVE, PUBLIC                             :: nb_trcobc    ! number of tracers with open BC 
     32   INTEGER  , SAVE, PUBLIC                             :: nb_trcsbc    ! number of tracers with surface BC 
     33   INTEGER  , SAVE, PUBLIC                             :: nb_trccbc    ! number of tracers with coastal BC 
    2934   INTEGER  , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: n_trc_indobc ! index of tracer with OBC data 
    3035   INTEGER  , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: n_trc_indsbc ! index of tracer with SBC data 
    3136   INTEGER  , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: n_trc_indcbc ! index of tracer with CBC data 
    32    INTEGER  , SAVE, PUBLIC                             :: ntra_obc     ! MAX( 1, nb_trcxxx ) to avoid compilation error with bounds checking 
    33    INTEGER  , SAVE, PUBLIC                             :: ntra_sbc     ! MAX( 1, nb_trcxxx ) to avoid compilation error with bounds checking 
    34    INTEGER  , SAVE, PUBLIC                             :: ntra_cbc     ! MAX( 1, nb_trcxxx ) to avoid compilation error with bounds checking 
    35    REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: rf_trofac   ! multiplicative factor for OBCtracer values 
    36    TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: sf_trcobc   ! structure of data input OBC (file informations, fields read) 
    37    REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: rf_trsfac   ! multiplicative factor for SBC tracer values 
    38    TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: sf_trcsbc   ! structure of data input SBC (file informations, fields read) 
    39    REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: rf_trcfac   ! multiplicative factor for CBC tracer values 
    40    TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: sf_trccbc   ! structure of data input CBC (file informations, fields read) 
     37   REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: rf_trsfac    ! multiplicative factor for SBC tracer values 
     38   TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: sf_trcsbc    ! structure of data input SBC (file informations, fields read) 
     39   REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: rf_trcfac    ! multiplicative factor for CBC tracer values 
     40   TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: sf_trccbc    ! structure of data input CBC (file informations, fields read) 
     41   REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:,:)  :: rf_trofac    ! multiplicative factor for OBCtracer values 
     42   TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:,:), TARGET  :: sf_trcobc    ! structure of data input OBC (file informations, fields read) 
     43   TYPE(MAP_POINTER), ALLOCATABLE, DIMENSION(:,:) :: nbmap_ptr   ! array of pointers to nbmap 
    4144 
    4245   !! * Substitutions 
    4346#  include "domzgr_substitute.h90" 
    4447   !!---------------------------------------------------------------------- 
    45    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     48   !! NEMO/OPA 3.6 , NEMO Consortium (2015) 
    4649   !! $Id$ 
    4750   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    6063      ! 
    6164      INTEGER,INTENT(IN) :: ntrc                           ! number of tracers 
    62       INTEGER            :: jl, jn                         ! dummy loop indices 
     65      INTEGER            :: jl, jn , ib, ibd, ii, ij, ik   ! dummy loop indices 
    6366      INTEGER            :: ierr0, ierr1, ierr2, ierr3     ! temporary integers 
    64       INTEGER            ::  ios                           ! Local integer output status for namelist read 
     67      INTEGER            :: ios                            ! Local integer output status for namelist read 
     68      INTEGER            :: nblen, igrd                    ! support arrays for BDY 
    6569      CHARACTER(len=100) :: clndta, clntrc 
    6670      ! 
    67       CHARACTER(len=100) :: cn_dir 
     71      CHARACTER(len=100) :: cn_dir_sbc, cn_dir_cbc, cn_dir_obc 
     72 
    6873      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i  ! local array of namelist informations on the fields to read 
    69       TYPE(FLD_N), DIMENSION(jpmaxtrc) :: sn_trcobc    ! open 
     74      TYPE(FLD_N), DIMENSION(jpmaxtrc,2) :: sn_trcobc  ! open 
     75      TYPE(FLD_N), DIMENSION(jpmaxtrc) :: sn_trcobc2   ! to read in multiple (2) open bdy 
    7076      TYPE(FLD_N), DIMENSION(jpmaxtrc) :: sn_trcsbc    ! surface 
    7177      TYPE(FLD_N), DIMENSION(jpmaxtrc) :: sn_trccbc    ! coastal 
    … …  
    7480      REAL(wp)   , DIMENSION(jpmaxtrc) :: rn_trcfac    ! multiplicative factor for tracer values 
    7581      !! 
    76       NAMELIST/namtrc_bc/ cn_dir, sn_trcobc, rn_trofac, sn_trcsbc, rn_trsfac, sn_trccbc, rn_trcfac  
     82      NAMELIST/namtrc_bc/ cn_dir_sbc, cn_dir_cbc, cn_dir_obc, sn_trcobc2, rn_trofac, sn_trcsbc, rn_trsfac, sn_trccbc, rn_trcfac 
     83#if defined key_bdy 
     84      NAMELIST/namtrc_bdy/ cn_trc_dflt, cn_trc, nn_trcdmp_bdy 
     85#endif 
    7786      !!---------------------------------------------------------------------- 
    7887      IF( nn_timing == 1 )  CALL timing_start('trc_bc_init') 
    7988      ! 
     89      IF( lwp ) THEN 
     90         WRITE(numout,*) ' ' 
     91         WRITE(numout,*) 'trc_bc_init : Tracers Boundary Conditions (BC)' 
     92         WRITE(numout,*) '~~~~~~~~~~~ ' 
     93      ENDIF 
    8094      !  Initialisation and local array allocation 
    8195      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0   
    … …  
    107121      n_trc_indcbc(:) = 0 
    108122      ! 
    109       DO jn = 1, ntrc 
    110          IF( ln_trc_obc(jn) ) THEN 
    111              nb_trcobc       = nb_trcobc + 1  
    112              n_trc_indobc(jn) = nb_trcobc  
    113          ENDIF 
    114          IF( ln_trc_sbc(jn) ) THEN 
    115              nb_trcsbc       = nb_trcsbc + 1 
    116              n_trc_indsbc(jn) = nb_trcsbc 
    117          ENDIF 
    118          IF( ln_trc_cbc(jn) ) THEN 
    119              nb_trccbc       = nb_trccbc + 1 
    120              n_trc_indcbc(jn) = nb_trccbc 
    121          ENDIF 
    122       ENDDO 
    123       ntra_obc = MAX( 1, nb_trcobc )   ! To avoid compilation error with bounds checking 
    124       IF( lwp ) WRITE(numout,*) ' ' 
    125       IF( lwp ) WRITE(numout,*) ' Number of passive tracers to be initialized with open boundary data :', nb_trcobc 
    126       IF( lwp ) WRITE(numout,*) ' ' 
    127       ntra_sbc = MAX( 1, nb_trcsbc )   ! To avoid compilation error with bounds checking 
    128       IF( lwp ) WRITE(numout,*) ' ' 
    129       IF( lwp ) WRITE(numout,*) ' Number of passive tracers to be initialized with surface boundary data :', nb_trcsbc 
    130       IF( lwp ) WRITE(numout,*) ' ' 
    131       ntra_cbc = MAX( 1, nb_trccbc )   ! To avoid compilation error with bounds checking 
    132       IF( lwp ) WRITE(numout,*) ' ' 
    133       IF( lwp ) WRITE(numout,*) ' Number of passive tracers to be initialized with coastal boundary data :', nb_trccbc 
    134       IF( lwp ) WRITE(numout,*) ' ' 
    135  
     123      ! Read Boundary Conditions Namelists 
    136124      REWIND( numnat_ref )              ! Namelist namtrc_bc in reference namelist : Passive tracer data structure 
    137125      READ  ( numnat_ref, namtrc_bc, IOSTAT = ios, ERR = 901) 
    … …  
    139127 
    140128      REWIND( numnat_cfg )              ! Namelist namtrc_bc in configuration namelist : Passive tracer data structure 
    141       READ  ( numnat_cfg, namtrc_bc, IOSTAT = ios, ERR = 902 ) 
    142 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_bc in configuration namelist', lwp ) 
    143       IF(lwm) WRITE ( numont, namtrc_bc ) 
    144  
    145       ! print some information for each  
     129#if defined key_bdy 
     130      DO ib = 1, nb_bdy 
     131#endif 
     132        READ  ( numnat_cfg, namtrc_bc, IOSTAT = ios, ERR = 902 ) 
     133902     IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_bc in configuration namelist', lwp ) 
     134        IF(lwm) WRITE ( numont, namtrc_bc ) 
     135#if defined key_bdy 
     136        sn_trcobc(:,ib)=sn_trcobc2(:) 
     137      ENDDO 
     138#endif 
     139 
     140#if defined key_bdy 
     141      REWIND( numnat_ref )              ! Namelist namtrc_bc in reference namelist : Passive tracer data structure 
     142      READ  ( numnat_ref, namtrc_bdy, IOSTAT = ios, ERR = 903) 
     143903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_bdy in reference namelist', lwp ) 
     144 
     145      REWIND( numnat_cfg )              ! Namelist namtrc_bc in configuration namelist : Passive tracer data structure 
     146      READ  ( numnat_cfg, namtrc_bdy, IOSTAT = ios, ERR = 904 ) 
     147904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_bdy in configuration namelist', lwp ) 
     148      IF(lwm) WRITE ( numont, namtrc_bdy ) 
     149      ! setup up preliminary informations for BDY structure 
     150      DO jn = 1, ntrc 
     151         DO ib = 1, nb_bdy 
     152            ! Set type of obc in BDY data structure (around here we may plug user override of obc type from nml) 
     153            IF ( ln_trc_obc(jn) ) THEN 
     154               trcdta_bdy(jn,ib)%cn_obc = TRIM( cn_trc(ib) ) 
     155            ELSE 
     156               trcdta_bdy(jn,ib)%cn_obc = TRIM( cn_trc_dflt(ib) ) 
     157            ENDIF 
     158            ! set damping use in BDY data structure 
     159            trcdta_bdy(jn,ib)%dmp = .false. 
     160            IF(nn_trcdmp_bdy(ib) .EQ. 1 .AND. ln_trc_obc(jn) ) trcdta_bdy(jn,ib)%dmp = .true. 
     161            IF(nn_trcdmp_bdy(ib) .EQ. 2 ) trcdta_bdy(jn,ib)%dmp = .true. 
     162            IF(trcdta_bdy(jn,ib)%cn_obc == 'frs' .AND. nn_trcdmp_bdy(ib) .NE. 0 )  & 
     163                & CALL ctl_stop( 'Use FRS OR relaxation' ) 
     164            IF (nn_trcdmp_bdy(ib) .LT. 0 .OR. nn_trcdmp_bdy(ib) .GT. 2)            & 
     165                & CALL ctl_stop( 'Not a valid option for nn_trcdmp_bdy. Allowed: 0,1,2.' ) 
     166         ENDDO 
     167      ENDDO 
     168 
     169#else 
     170      ! Force all tracers OBC to false if bdy not used 
     171      ln_trc_obc = .false. 
     172#endif 
     173      ! compose BC data indexes 
     174      DO jn = 1, ntrc 
     175         IF( ln_trc_obc(jn) ) THEN 
     176             nb_trcobc       = nb_trcobc + 1  ; n_trc_indobc(jn) = nb_trcobc 
     177         ENDIF 
     178         IF( ln_trc_sbc(jn) ) THEN 
     179             nb_trcsbc       = nb_trcsbc + 1  ; n_trc_indsbc(jn) = nb_trcsbc 
     180         ENDIF 
     181         IF( ln_trc_cbc(jn) ) THEN 
     182             nb_trccbc       = nb_trccbc + 1  ; n_trc_indcbc(jn) = nb_trccbc 
     183         ENDIF 
     184      ENDDO 
     185 
     186      ! Print summmary of Boundary Conditions 
    146187      IF( lwp ) THEN 
    147          DO jn = 1, ntrc 
    148             IF( ln_trc_obc(jn) )  THEN     
    149                clndta = TRIM( sn_trcobc(jn)%clvar )  
    150                IF(lwp) WRITE(numout,*) 'Preparing to read OBC data file for passive tracer number :', jn, ' name : ', clndta, &  
    151                &               ' multiplicative factor : ', rn_trofac(jn) 
    152             ENDIF 
    153             IF( ln_trc_sbc(jn) )  THEN     
    154                clndta = TRIM( sn_trcsbc(jn)%clvar )  
    155                IF(lwp) WRITE(numout,*) 'Preparing to read SBC data file for passive tracer number :', jn, ' name : ', clndta, &  
    156                &               ' multiplicative factor : ', rn_trsfac(jn) 
    157             ENDIF 
    158             IF( ln_trc_cbc(jn) )  THEN     
    159                clndta = TRIM( sn_trccbc(jn)%clvar )  
    160                IF(lwp) WRITE(numout,*) 'Preparing to read CBC data file for passive tracer number :', jn, ' name : ', clndta, &  
    161                &               ' multiplicative factor : ', rn_trcfac(jn) 
    162             ENDIF 
    163          END DO 
    164       ENDIF