Changeset 8059


Ignore:
Timestamp:
2017-05-23T10:32:39+02:00 (3 years ago)
Author:
jgraham
Message:

Merged branches required for AMM15 simulations, see ticket #1904.
Merged branches include:
branches/UKMO/CO6_KD490
branches/UKMO/CO6_Restartable_Tidal_Analysis
branches/UKMO/AMM15_v3_6_STABLE

Location:
branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC
Files:
5 added
32 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r8058 r8059  
    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 
  • branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r8058 r8059  
    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 
  • branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90

    r8058 r8059  
    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/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r8058 r8059  
    213213             dta_bdy(ib_bdy)%ll_u3d = .true. 
    214214             dta_bdy(ib_bdy)%ll_v3d = .true. 
     215          CASE('neumann')  
     216             IF(lwp) WRITE(numout,*) '      Neumann conditions'  
     217             dta_bdy(ib_bdy)%ll_u3d = .false.  
     218             dta_bdy(ib_bdy)%ll_v3d = .false.  
     219          CASE('zerograd')  
     220             IF(lwp) WRITE(numout,*) '      Zero gradient for baroclinic velocities'  
     221             dta_bdy(ib_bdy)%ll_u3d = .false.  
     222             dta_bdy(ib_bdy)%ll_v3d = .false. 
    215223          CASE('zero') 
    216224             IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
     
    10871095            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    10881096               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 
     1097               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) * 0.5 & 
     1098                                            &  *(10./ FLOAT(nn_rimwidth(ib_bdy))) ) ! JGraham:modified for rim=15 
     1099!               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation 
    10901100!               idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.  ! quadratic 
    10911101!               idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy))       ! linear 
  • branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90

    r8058 r8059  
    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/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90

    r8058 r8059  
    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/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90

    r8058 r8059  
    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/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r8058 r8059  
    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 
    4648   USE iom 
    4749   USE ioipsl 
    4850   USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities      
     51   USE eosbn2         ! equation of state                (eos_bn2 routine) 
     52 
    4953 
    5054#if defined key_lim2 
     
    132136      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace 
    133137      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace 
     138      REAL(wp), POINTER, DIMENSION(:,:,:) :: zrhd , zrhop  ! 3D workspace 
    134139      !!---------------------------------------------------------------------- 
    135140      !  
     
    138143      CALL wrk_alloc( jpi , jpj      , z2d ) 
    139144      CALL wrk_alloc( jpi , jpj, jpk , z3d ) 
     145      CALL wrk_alloc( jpi , jpj, jpk , zrhd , zrhop ) 
    140146      ! 
    141147      ! Output the initial state and forcings 
     
    376382         CALL iom_put( "v_salttr", 0.5 * z2d )            !  heat transport in j-direction 
    377383      ENDIF 
     384 
     385      IF( iom_use("rhop") ) THEN 
     386         CALL eos( tsn, zrhd, zrhop, fsdept_n(:,:,:) )       ! now in situ and potential density 
     387         zrhop(:,:,jpk) = 0._wp 
     388         CALL iom_put( 'rhop', zrhop ) 
     389      ENDIF 
     390 
    378391      ! 
    379392      CALL wrk_dealloc( jpi , jpj      , z2d ) 
    380393      CALL wrk_dealloc( jpi , jpj, jpk , z3d ) 
     394      CALL wrk_dealloc( jpi , jpj, jpk , zrhd , zrhop ) 
     395      ! 
     396      ! If we want tmb values  
     397 
     398      IF (ln_diatmb) THEN 
     399         CALL dia_tmb 
     400      ENDIF 
     401      IF (ln_dia25h) THEN 
     402         CALL dia_25h( kt ) 
     403      ENDIF 
    381404      ! 
    382405      IF( nn_timing == 1 )   CALL timing_stop('dia_wri') 
  • branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r8058 r8059  
    136136      USE ioipsl 
    137137      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,   & 
     138         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstdate, ln_rstart , nn_rstctl,   & 
    139139         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   & 
    140140         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler 
     
    174174         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir 
    175175         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart 
     176         WRITE(numout,*) '      datestamping of restarts        ln_rstdate  = ', ln_rstdate 
    176177         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler 
    177178         WRITE(numout,*) '      control of time step            nn_rstctl  = ', nn_rstctl 
  • branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r8058 r8059  
    3131   USE wrk_nemo        ! Memory allocation 
    3232   USE timing          ! Timing 
     33   USE iom    ! slwa 
    3334 
    3435   IMPLICIT NONE 
  • branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90

    r8058 r8059  
    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/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90

    r8058 r8059  
    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      ! 
     
    134170      END SELECT 
    135171      ! 
     172 
     173#if defined key_bdy  
     174      ! restore velocity masks at points outside boundary  
     175      un(:,:,:) = un(:,:,:) * umask(:,:,:)  
     176      vn(:,:,:) = vn(:,:,:) * vmask(:,:,:)  
     177#endif   
     178 
    136179      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==! 
    137180         DO jj = 2, jpjm1 
  • branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r8058 r8059  
    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 
     
    144145      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices 
    145146      INTEGER  ::   ikbu, ikbv, noffset      ! local integers 
     147      REAL(wp) ::   zmdi 
    146148      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars 
    147149      REAL(wp) ::   zx1, zy1, zx2, zy2         !   -      - 
     
    169171      CALL wrk_alloc( jpi, jpj, zhf ) 
    170172      ! 
     173      zmdi=1.e+20                               !  missing data indicator for masking 
    171174      !                                         !* Local constant initialization 
    172175      z1_12 = 1._wp / 12._wp  
     
    926929      CALL wrk_dealloc( jpi, jpj, zhf ) 
    927930      ! 
     931      IF ( ln_diatmb ) THEN 
     932         CALL iom_put( "baro_u" , un_b*umask(:,:,1)+zmdi*(1-umask(:,:,1 ) ) )  ! Barotropic  U Velocity 
     933         CALL iom_put( "baro_v" , vn_b*vmask(:,:,1)+zmdi*(1-vmask(:,:,1 ) ) )  ! Barotropic  V Velocity 
     934      ENDIF 
    928935      IF( nn_timing == 1 )  CALL timing_stop('dyn_spg_ts') 
    929936      ! 
  • branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/ICB/icbrst.F90

    r8058 r8059  
    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/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90

    r8058 r8059  
    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 
  • branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r8058 r8059  
    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/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r8058 r8059  
    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/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90

    r8058 r8059  
    2828   PUBLIC sbc_flx       ! routine called by step.F90 
    2929 
    30    INTEGER , PARAMETER ::   jpfld   = 5   ! maximum number of files to read  
     30   INTEGER , PARAMETER ::   jpfld   = 6   ! maximum number of files to read  
    3131   INTEGER , PARAMETER ::   jp_utau = 1   ! index of wind stress (i-component) file 
    3232   INTEGER , PARAMETER ::   jp_vtau = 2   ! index of wind stress (j-component) file 
     
    3434   INTEGER , PARAMETER ::   jp_qsr  = 4   ! index of solar heat file 
    3535   INTEGER , PARAMETER ::   jp_emp  = 5   ! index of evaporation-precipation file 
     36   INTEGER , PARAMETER ::   jp_press = 6  ! index of pressure for UKMO shelf fluxes 
    3637   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read) 
     38   LOGICAL , PUBLIC    ::   ln_shelf_flx = .FALSE. ! UKMO SHELF specific flux flag 
     39   INTEGER             ::   jpfld_local   ! maximum number of files to read (locally modified depending on ln_shelf_flx)  
    3740 
    3841   !! * Substitutions 
     
    8285      REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient 
    8386      REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables 
     87      REAL     ::   cs                    ! UKMO SHELF: Friction co-efficient at surface 
     88      REAL     ::   totwindspd            ! UKMO SHELF: Magnitude of wind speed vector 
     89 
     90      REAL(wp) ::   rhoa  = 1.22         ! Air density kg/m3 
     91      REAL(wp) ::   cdrag = 1.5e-3       ! drag coefficient  
    8492      !! 
    8593      CHARACTER(len=100) ::  cn_dir                               ! Root directory for location of flx files 
    8694      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 
     95      TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp, sn_press  !  informations about the fields to be read 
     96      LOGICAL     ::   ln_foam_flx  = .FALSE.                     ! UKMO FOAM specific flux flag 
     97      NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp,   & 
     98      &                    ln_foam_flx, sn_press, ln_shelf_flx 
    8999      !!--------------------------------------------------------------------- 
    90100      ! 
     
    109119         slf_i(jp_emp ) = sn_emp 
    110120         ! 
    111          ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure 
     121            ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure 
     122            IF( ln_shelf_flx ) slf_i(jp_press) = sn_press 
     123    
     124            ! define local jpfld depending on shelf_flx logical 
     125            IF( ln_shelf_flx ) THEN 
     126               jpfld_local = jpfld 
     127            ELSE 
     128               jpfld_local = jpfld-1 
     129            ENDIF 
     130            ! 
    112131         IF( ierror > 0 ) THEN    
    113132            CALL ctl_stop( 'sbc_flx: unable to allocate sf structure' )   ;   RETURN   
    114133         ENDIF 
    115          DO ji= 1, jpfld 
     134         DO ji= 1, jpfld_local 
    116135            ALLOCATE( sf(ji)%fnow(jpi,jpj,1) ) 
    117136            IF( slf_i(ji)%ln_tint ) ALLOCATE( sf(ji)%fdta(jpi,jpj,1,2) ) 
     
    132151         ENDIF 
    133152!CDIR COLLAPSE 
     153            !!UKMO SHELF effect of atmospheric pressure on SSH 
     154            ! If using ln_apr_dyn, this is done there so don't repeat here. 
     155            IF( ln_shelf_flx .AND. .NOT. ln_apr_dyn) THEN 
     156               DO jj = 1, jpjm1 
     157                  DO ji = 1, jpim1 
     158                     apgu(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji+1,jj,1)-sf(jp_press)%fnow(ji,jj,1))/e1u(ji,jj) 
     159                     apgv(ji,jj) = (-1.0/rau0)*(sf(jp_press)%fnow(ji,jj+1,1)-sf(jp_press)%fnow(ji,jj,1))/e2v(ji,jj) 
     160                  END DO 
     161               END DO 
     162            ENDIF ! ln_shelf_flx 
     163       
    134164         DO jj = 1, jpj                                           ! set the ocean fluxes from read fields 
    135165            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) 
     166                   IF( ln_shelf_flx ) THEN 
     167                      !! UKMO SHELF - need atmospheric pressure to calculate Haney forcing 
     168                      pressnow(ji,jj) = sf(jp_press)%fnow(ji,jj,1) 
     169                      !! UKMO SHELF flux files contain wind speed not wind stress 
     170                      totwindspd = sqrt((sf(jp_utau)%fnow(ji,jj,1))**2.0 + (sf(jp_vtau)%fnow(ji,jj,1))**2.0) 
     171                      cs = 0.63 + (0.066 * totwindspd) 
     172                      utau(ji,jj) = cs * (rhoa/rau0) * sf(jp_utau)%fnow(ji,jj,1) * totwindspd 
     173                      vtau(ji,jj) = cs * (rhoa/rau0) * sf(jp_vtau)%fnow(ji,jj,1) * totwindspd 
     174                   ELSE 
     175                      utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) 
     176                      vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) 
     177                   ENDIF 
     178                   qsr (ji,jj) = sf(jp_qsr )%fnow(ji,jj,1) 
     179                   IF( ln_foam_flx .OR. ln_shelf_flx ) THEN 
     180                      !! UKMO FOAM flux files contain non-solar heat flux (qns) rather than total heat flux (qtot) 
     181                      qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) 
     182                      !! UKMO FOAM flux files contain the net DOWNWARD freshwater flux P-E rather then E-P 
     183                      emp (ji,jj) = -1. * sf(jp_emp )%fnow(ji,jj,1) 
     184                   ELSE 
     185                      qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) 
     186                      emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1) 
     187                   ENDIF 
    140188            END DO 
    141189         END DO 
     
    143191         qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST 
    144192         ! 
     193    
     194            !! UKMO FOAM wind fluxes need lbc_lnk calls owing to a bug in interp.exe 
     195            IF( ln_foam_flx ) THEN 
     196               CALL lbc_lnk( utau(:,:), 'U', -1. ) 
     197               CALL lbc_lnk( vtau(:,:), 'V', -1. ) 
     198            ENDIF 
     199     
    145200         !                                                        ! module of wind stress and wind speed at T-point 
    146201         zcoef = 1. / ( zrhoa * zcdrag ) 
     
    162217            WRITE(numout,*)  
    163218            WRITE(numout,*) '        read daily momentum, heat and freshwater fluxes OK' 
    164             DO jf = 1, jpfld 
     219            DO jf = 1, jpfld_local 
    165220               IF( jf == jp_utau .OR. jf == jp_vtau )   zfact =     1. 
    166221               IF( jf == jp_qtot .OR. jf == jp_qsr  )   zfact =     0.1 
  • branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90

    r8058 r8059  
    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/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/tide.h90

    r4292 r8059  
    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/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/SBC/tide_mod.F90

    r8058 r8059  
    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/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r8058 r8059  
    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/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90

    r8058 r8059  
    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/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

    r8058 r8059  
    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      !!---------------------------------------------------------------------- 
     
    120121      IF( lk_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries 
    121122#endif 
    122   
     123 
     124#if ( ! defined key_lim3 && ! defined key_lim2 && ! key_cice ) 
     125      IF ( kt == nit000 ) warn_1=0 
     126      warn_2=0 
     127      DO jk = 1, jpkm1 
     128         DO jj = 1, jpj 
     129            DO ji = 1, jpi 
     130               IF ( tsa(ji,jj,jk,jp_tem) .lt. 0.0 ) THEN 
     131                  ! calculate freezing point 
     132                  zfreeze = ( -0.0575_wp + 1.710523E-3 * Sqrt(Abs(tsn(ji,jj,jk,jp_sal)))   &  
     133                            - 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) ) 
     134                  IF ( tsa(ji,jj,jk,jp_tem) .lt. zfreeze ) THEN 
     135                     tsa(ji,jj,jk,jp_tem)=zfreeze 
     136                     warn_2=1 
     137                  ENDIF 
     138               ENDIF 
     139            END DO 
     140         END DO 
     141      END DO 
     142      CALL mpp_max(warn_1) 
     143      CALL mpp_max(warn_2) 
     144      IF ( (warn_1 == 0) .and. (warn_2 /= 0) ) THEN 
     145         IF(lwp) THEN 
     146            CALL ctl_warn( ' Temperatures dropping below freezing point, ', & 
     147                      &    ' being forced to freezing point, no longer conservative' )  
     148         ENDIF 
     149         warn_1=1 
     150      ENDIF 
     151#endif 
     152 
    123153      ! set time step size (Euler/Leapfrog) 
    124154      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dtra(:) =     rdttra(:)      ! at nit000             (Euler) 
  • branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r8058 r8059  
    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/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r8058 r8059  
    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 
     
    138139 
    139140!!gm      IF( .NOT.ln_traqsr )   qsr(:,:) = 0.e0   ! no solar radiation penetration 
    140       IF( .NOT.ln_traqsr ) THEN     ! no solar radiation penetration 
     141      IF( .NOT.ln_traqsr .and. .NOT.ln_tradwl ) THEN     ! no solar radiation penetration 
    141142         qns(:,:) = qns(:,:) + qsr(:,:)      ! total heat flux in qns 
    142143         qsr(:,:) = 0.e0                     ! qsr set to zero 
  • branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90

    r8058 r8059  
    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/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

    r8058 r8059  
    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 
     
    2728 
    2829   PUBLIC   zdf_mxl       ! called by step.F90 
    29    PUBLIC   zdf_mxl_alloc ! Used in zdf_tke_init 
    3030 
    3131   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln    !: number of level in the mixed layer (used by TOP) 
     
    3333   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m] 
    3434   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m] 
     35   REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:) ::   hmld_zint  !: vertically-interpolated mixed layer depth   [m]  
     36   REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:) ::   htc_mld    ! Heat content of hmld_zint 
     37   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    :: ll_found   ! Is T_b to be found by interpolation ?  
     38   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: ll_belowml ! Flag points below mixed layer when ll_found=F 
    3539 
    3640   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
    3741   REAL(wp)         ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth 
     42 
     43   TYPE, PUBLIC :: MXL_ZINT   !: Structure for MLD defs 
     44      INTEGER   :: mld_type   ! mixed layer type      
     45      REAL(wp)  :: zref       ! depth of initial T_ref 
     46      REAL(wp)  :: dT_crit    ! Critical temp diff 
     47      REAL(wp)  :: iso_frac   ! Fraction of rn_dT_crit used 
     48   END TYPE MXL_ZINT 
    3849 
    3950   !! * Substitutions 
     
    5263      zdf_mxl_alloc = 0      ! set to zero if no array to be allocated 
    5364      IF( .NOT. ALLOCATED( nmln ) ) THEN 
    54          ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc ) 
     65         ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj),       & 
     66        &          htc_mld(jpi,jpj),                                                                    & 
     67        &          ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc ) 
    5568         ! 
    5669         IF( lk_mpp             )   CALL mpp_sum ( zdf_mxl_alloc ) 
     
    128141            iikn = nmln(ji,jj) 
    129142            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 
     143            hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,imkt )            )  * ssmask(ji,jj)    ! Turbocline depth  
     144            hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,MAX( imkt,nla10 ) ) ) * ssmask(ji,jj)    ! Mixed layer depth 
     145            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 
    133146         END DO 
    134147      END DO 
     
    138151      ENDIF 
    139152       
     153      ! Vertically-interpolated mixed-layer depth diagnostic 
     154      CALL zdf_mxl_zint( kt ) 
     155 
    140156      IF(ln_ctl)   CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 ) 
    141157      ! 
     
    145161      ! 
    146162   END SUBROUTINE zdf_mxl 
     163 
     164   SUBROUTINE zdf_mxl_zint_mld( sf )  
     165      !!----------------------------------------------------------------------------------  
     166      !!                    ***  ROUTINE zdf_mxl_zint_mld  ***  
     167      !                                                                         
     168      !   Calculate vertically-interpolated mixed layer depth diagnostic.  
     169      !             
     170      !   This routine can calculate the mixed layer depth diagnostic suggested by 
     171      !   Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate 
     172      !   vertically-interpolated mixed-layer depth diagnostics with other parameter 
     173      !   settings set in the namzdf_mldzint namelist.   
     174      !  
     175      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the   
     176      !   density has increased by an amount equivalent to a temperature difference of   
     177      !   0.8C at the surface.  
     178      !  
     179      !   For other values of mld_type the mixed layer is calculated as the depth at   
     180      !   which the temperature differs by 0.8C from the surface temperature.   
     181      !                                                                         
     182      !   David Acreman, Daley Calvert                                       
     183      !  
     184      !!-----------------------------------------------------------------------------------  
     185 
     186      TYPE(MXL_ZINT), INTENT(in)  :: sf 
     187 
     188      ! Diagnostic criteria 
     189      INTEGER   :: nn_mld_type   ! mixed layer type      
     190      REAL(wp)  :: rn_zref       ! depth of initial T_ref 
     191      REAL(wp)  :: rn_dT_crit    ! Critical temp diff 
     192      REAL(wp)  :: rn_iso_frac   ! Fraction of rn_dT_crit used 
     193 
     194      ! Local variables 
     195      REAL(wp), PARAMETER :: zepsilon = 1.e-30          ! local small value 
     196      INTEGER, POINTER, DIMENSION(:,:) :: ikmt          ! number of active tracer levels  
     197      INTEGER, POINTER, DIMENSION(:,:) :: ik_ref        ! index of reference level  
     198      INTEGER, POINTER, DIMENSION(:,:) :: ik_iso        ! index of last uniform temp level  
     199      REAL, POINTER, DIMENSION(:,:,:)  :: zT            ! Temperature or density  
     200      REAL, POINTER, DIMENSION(:,:)    :: ppzdep        ! depth for use in calculating d(rho)  
     201      REAL, POINTER, DIMENSION(:,:)    :: zT_ref        ! reference temperature  
     202      REAL    :: zT_b                                   ! base temperature  
     203      REAL, POINTER, DIMENSION(:,:,:)  :: zdTdz         ! gradient of zT  
     204      REAL, POINTER, DIMENSION(:,:,:)  :: zmoddT        ! Absolute temperature difference  
     205      REAL    :: zdz                                    ! depth difference  
     206      REAL    :: zdT                                    ! temperature difference  
     207      REAL, POINTER, DIMENSION(:,:)    :: zdelta_T      ! difference critereon  
     208      REAL, POINTER, DIMENSION(:,:)    :: zRHO1, zRHO2  ! Densities  
     209      INTEGER :: ji, jj, jk                             ! loop counter  
     210 
     211      !!-------------------------------------------------------------------------------------  
     212      !   
     213      CALL wrk_alloc( jpi, jpj, ikmt, ik_ref, ik_iso)  
     214      CALL wrk_alloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 )  
     215      CALL wrk_alloc( jpi, jpj, jpk, zT, zdTdz, zmoddT )  
     216 
     217      ! Unpack structure 
     218      nn_mld_type = sf%mld_type 
     219      rn_zref     = sf%zref 
     220      rn_dT_crit  = sf%dT_crit 
     221      rn_iso_frac = sf%iso_frac 
     222 
     223      ! Set the mixed layer depth criterion at each grid point  
     224      IF( nn_mld_type == 0 ) THEN 
     225         zdelta_T(:,:) = rn_dT_crit 
     226         zT(:,:,:) = rhop(:,:,:) 
     227      ELSE IF( nn_mld_type == 1 ) THEN 
     228         ppzdep(:,:)=0.0  
     229         call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) )  
     230! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST  
     231! [assumes number of tracers less than number of vertical levels]  
     232         zT(:,:,1:jpts)=tsn(:,:,1,1:jpts)  
     233         zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit  
     234         CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) )  
     235         zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0  
     236         ! RHO from eos (2d version) doesn't calculate north or east halo:  
     237         CALL lbc_lnk( zdelta_T, 'T', 1. )  
     238         zT(:,:,:) = rhop(:,:,:)  
     239      ELSE  
     240         zdelta_T(:,:) = rn_dT_crit                       
     241         zT(:,:,:) = tsn(:,:,:,jp_tem)                            
     242      END IF  
     243 
     244      ! Calculate the gradient of zT and absolute difference for use later  
     245      DO jk = 1 ,jpk-2  
     246         zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / fse3w(:,:,jk+1)  
     247         zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )  
     248      END DO  
     249 
     250      ! Find density/temperature at the reference level (Kara et al use 10m).           
     251      ! ik_ref is the index of the box centre immediately above or at the reference level  
     252      ! Find rn_zref in the array of model level depths and find the ref     
     253      ! density/temperature by linear interpolation.                                    
     254      DO jk = jpkm1, 2, -1  
     255         WHERE ( fsdept(:,:,jk) > rn_zref )  
     256           ik_ref(:,:) = jk - 1  
     257           zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - fsdept(:,:,jk-1) )  
     258         END WHERE  
     259      END DO  
     260 
     261      ! If the first grid box centre is below the reference level then use the  
     262      ! top model level to get zT_ref  
     263      WHERE ( fsdept(:,:,1) > rn_zref )   
     264         zT_ref = zT(:,:,1)  
     265         ik_ref = 1  
     266      END WHERE  
     267 
     268      ! The number of active tracer levels is 1 less than the number of active w levels  
     269      ikmt(:,:) = mbathy(:,:) - 1  
     270 
     271      ! Initialize / reset 
     272      ll_found(:,:) = .false. 
     273 
     274      IF ( rn_iso_frac - zepsilon > 0. ) THEN 
     275         ! Search for a uniform density/temperature region where adjacent levels           
     276         ! differ by less than rn_iso_frac * deltaT.                                       
     277         ! ik_iso is the index of the last level in the uniform layer   
     278         ! ll_found indicates whether the mixed layer depth can be found by interpolation  
     279         ik_iso(:,:)   = ik_ref(:,:)  
     280         DO jj = 1, nlcj  
     281            DO ji = 1, nlci  
     282!CDIR NOVECTOR  
     283               DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1  
     284                  IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN  
     285                     ik_iso(ji,jj)   = jk  
     286                     ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )  
     287                     EXIT  
     288                  END IF  
     289               END DO  
     290            END DO  
     291         END DO  
     292 
     293         ! Use linear interpolation to find depth of mixed layer base where possible  
     294         hmld_zint(:,:) = rn_zref  
     295         DO jj = 1, jpj  
     296            DO ji = 1, jpi  
     297               IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN  
     298                  zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )  
     299                  hmld_zint(ji,jj) = fsdept(ji,jj,ik_iso(ji,jj)) + zdz  
     300               END IF  
     301            END DO  
     302         END DO  
     303      END IF 
     304 
     305      ! If ll_found = .false. then calculate MLD using difference of zdelta_T     
     306      ! from the reference density/temperature  
     307  
     308! Prevent this section from working on land points  
     309      WHERE ( tmask(:,:,1) /= 1.0 )  
     310         ll_found = .true.  
     311      END WHERE  
     312  
     313      DO jk=1, jpk  
     314         ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:)   
     315      END DO  
     316  
     317! Set default value where interpolation cannot be used (ll_found=false)   
     318      DO jj = 1, jpj  
     319         DO ji = 1, jpi  
     320            IF ( .not. ll_found(ji,jj) )  hmld_zint(ji,jj) = fsdept(ji,jj,ikmt(ji,jj))  
     321         END DO  
     322      END DO  
     323 
     324      DO jj = 1, jpj  
     325         DO ji = 1, jpi  
     326!CDIR NOVECTOR  
     327            DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj)  
     328               IF ( ll_found(ji,jj) ) EXIT  
     329               IF ( ll_belowml(ji,jj,jk) ) THEN                 
     330                  zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) )  
     331                  zdT  = zT_b - zT(ji,jj,jk-1)                                       
     332                  zdz  = zdT / zdTdz(ji,jj,jk-1)                                        
     333                  hmld_zint(ji,jj) = fsdept(ji,jj,jk-1) + zdz  
     334                  EXIT                                                    
     335               END IF  
     336            END DO  
     337         END DO  
     338      END DO  
     339 
     340      hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1)  
     341      !   
     342      CALL wrk_dealloc( jpi, jpj, ikmt, ik_ref, ik_iso)  
     343      CALL wrk_dealloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 )  
     344      CALL wrk_dealloc( jpi,jpj, jpk, zT, zdTdz, zmoddT )  
     345      !  
     346   END SUBROUTINE zdf_mxl_zint_mld 
     347 
     348   SUBROUTINE zdf_mxl_zint_htc( kt ) 
     349      !!---------------------------------------------------------------------- 
     350      !!                  ***  ROUTINE zdf_mxl_zint_htc  *** 
     351      !!  
     352      !! ** Purpose :    
     353      !! 
     354      !! ** Method  :    
     355      !!---------------------------------------------------------------------- 
     356 
     357      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     358 
     359      INTEGER :: ji, jj, jk 
     360      INTEGER :: ikmax 
     361      REAL(wp) :: zc, zcoef 
     362      ! 
     363      INTEGER,  ALLOCATABLE, DIMENSION(:,:) ::   ilevel 
     364      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zthick_0, zthick 
     365 
     366      !!---------------------------------------------------------------------- 
     367 
     368      IF( .NOT. ALLOCATED(ilevel) ) THEN 
     369         ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), & 
     370         &         zthick(jpi,jpj), STAT=ji ) 
     371         IF( lk_mpp  )   CALL mpp_sum(ji) 
     372         IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' ) 
     373      ENDIF 
     374 
     375      ! Find last whole model T level above the MLD 
     376      ilevel(:,:)   = 0 
     377      zthick_0(:,:) = 0._wp 
     378 
     379      DO jk = 1, jpkm1   
     380         DO jj = 1, jpj 
     381            DO ji = 1, jpi                     
     382               zthick_0(ji,jj) = zthick_0(ji,jj) + fse3t(ji,jj,jk) 
     383               IF( zthick_0(ji,jj) < hmld_zint(ji,jj) )   ilevel(ji,jj) = jk 
     384            END DO 
     385         END DO 
     386         WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2) 
     387         WRITE(numout,*) 'fsdepw(jk+1 =',jk+1,') =',fsdepw(2,2,jk+1) 
     388      END DO 
     389 
     390      ! Surface boundary condition 
     391      IF( lk_vvl ) THEN   ;   zthick(:,:) = 0._wp       ;   htc_mld(:,:) = 0._wp                                    
     392      ELSE                ;   zthick(:,:) = sshn(:,:)   ;   htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)    
     393      ENDIF 
     394 
     395      ! Deepest whole T level above the MLD 
     396      ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 ) 
     397 
     398      ! Integration down to last whole model T level 
     399      DO jk = 1, ikmax 
     400         DO jj = 1, jpj 
     401            DO ji = 1, jpi 
     402               zc = fse3t(ji,jj,jk) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1  )  )    ! 0 below ilevel 
     403               zthick(ji,jj) = zthick(ji,jj) + zc 
     404               htc_mld(ji,jj) = htc_mld(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 
     405            END DO 
     406         END DO 
     407      END DO 
     408 
     409      ! Subsequent partial T level 
     410      zthick(:,:) = hmld_zint(:,:) - zthick(:,:)   !   remaining thickness to reach MLD 
     411 
     412      DO jj = 1, jpj 
     413         DO ji = 1, jpi 
     414            htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem)  &  
     415      &                      * MIN( fse3t(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1) 
     416         END DO 
     417      END DO 
     418 
     419      WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2) 
     420 
     421      ! Convert to heat content 
     422      zcoef = rau0 * rcp 
     423      htc_mld(:,:) = zcoef * htc_mld(:,:) 
     424 
     425   END SUBROUTINE zdf_mxl_zint_htc 
     426 
     427   SUBROUTINE zdf_mxl_zint( kt ) 
     428      !!---------------------------------------------------------------------- 
     429      !!                  ***  ROUTINE zdf_mxl_zint  *** 
     430      !!  
     431      !! ** Purpose :    
     432      !! 
     433      !! ** Method  :    
     434      !!---------------------------------------------------------------------- 
     435 
     436      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     437 
     438      INTEGER :: ios 
     439      INTEGER :: jn 
     440 
     441      INTEGER :: nn_mld_diag = 0    ! number of diagnostics 
     442 
     443      CHARACTER(len=1) :: cmld 
     444 
     445      TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 
     446      TYPE(MXL_ZINT), SAVE, DIMENSION(5) ::   mld_diags 
     447 
     448      NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 
     449 
     450      !!---------------------------------------------------------------------- 
     451       
     452      IF( kt == nit000 ) THEN 
     453         REWIND( numnam_ref )              ! Namelist namzdf_mldzint in reference namelist  
     454         READ  ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901) 
     455901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist', lwp ) 
     456 
     457         REWIND( numnam_cfg )              ! Namelist namzdf_mldzint in configuration namelist  
     458         READ  ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 ) 
     459902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist', lwp ) 
     460         IF(lwm) WRITE ( numond, namzdf_mldzint ) 
     461 
     462         IF( nn_mld_diag > 5 )   CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' ) 
     463 
     464         mld_diags(1) = sn_mld1 
     465         mld_diags(2) = sn_mld2 
     466         mld_diags(3) = sn_mld3 
     467         mld_diags(4) = sn_mld4 
     468         mld_diags(5) = sn_mld5 
     469 
     470         IF( nn_mld_diag > 0 ) THEN 
     471            WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================' 
     472            WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)' 
     473            DO jn = 1, nn_mld_diag 
     474               WRITE(numout,*) 'MLD criterion',jn,':' 
     475               WRITE(numout,*) '    nn_mld_type =', mld_diags(jn)%mld_type 
     476               WRITE(numout,*) '    rn_zref ='    , mld_diags(jn)%zref 
     477               WRITE(numout,*) '    rn_dT_crit =' , mld_diags(jn)%dT_crit 
     478               WRITE(numout,*) '    rn_iso_frac =', mld_diags(jn)%iso_frac 
     479            END DO 
     480            WRITE(numout,*) '====================================================================' 
     481         ENDIF 
     482      ENDIF 
     483 
     484      IF( nn_mld_diag > 0 ) THEN 
     485         DO jn = 1, nn_mld_diag 
     486            WRITE(cmld,'(I1)') jn 
     487            IF( iom_use( "mldzint_"//cmld ) .OR. iom_use( "mldhtc_"//cmld ) ) THEN 
     488               CALL zdf_mxl_zint_mld( mld_diags(jn) ) 
     489 
     490               IF( iom_use( "mldzint_"//cmld ) ) THEN 
     491                  CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) ) 
     492               ENDIF 
     493 
     494               IF( iom_use( "mldhtc_"//cmld ) )  THEN 
     495                  CALL zdf_mxl_zint_htc( kt ) 
     496                  CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:)   ) 
     497               ENDIF 
     498            ENDIF 
     499         END DO 
     500      ENDIF 
     501 
     502   END SUBROUTINE zdf_mxl_zint 
    147503 
    148504   !!====================================================================== 
  • branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r8058 r8059  
    8585   USE stopar 
    8686   USE stopts 
     87   USE diatmb          ! Top,middle,bottom output 
     88   USE dia25h          ! 25h mean output 
    8789 
    8890   IMPLICIT NONE 
     
    475477      IF( lk_asminc     )   CALL asm_inc_init   ! Initialize assimilation increments 
    476478      IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler 
     479                            CALL dia_tmb_init  ! TMB outputs 
     480                            CALL dia_25h_init  ! 25h mean  outputs 
    477481      ! 
    478482   END SUBROUTINE nemo_init 
     
    630634      USE ldftra_oce, ONLY: ldftra_oce_alloc 
    631635      USE trc_oce   , ONLY: trc_oce_alloc 
     636      USE diainsitutem, ONLY: insitu_tem_alloc 
    632637#if defined key_diadct  
    633638      USE diadct    , ONLY: diadct_alloc  
     
    646651      ierr = ierr + ldftra_oce_alloc()          ! ocean lateral  physics : tracers 
    647652      ierr = ierr + zdf_oce_alloc   ()          ! ocean vertical physics 
     653      ierr = ierr + insitu_tem_alloc() 
    648654      ! 
    649655      ierr = ierr + trc_oce_alloc   ()          ! shared TRC / TRA arrays 
  • branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/step.F90

    r8058 r8059  
    255255                             CALL tra_sbc    ( kstp )       ! surface boundary condition 
    256256      IF( ln_traqsr      )   CALL tra_qsr    ( kstp )       ! penetrative solar radiation qsr 
     257      IF( ln_tradwl      )   CALL tra_dwl    ( kstp )       ! Polcoms Style Short Wave Radiation  
    257258      IF( ln_trabbc      )   CALL tra_bbc    ( kstp )       ! bottom heat flux 
    258259      IF( lk_trabbl      )   CALL tra_bbl    ( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme 
     
    337338      IF( lk_vvl           )   CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors 
    338339      ! 
    339       IF( lrst_oce         )   CALL rst_write( kstp )       ! write output ocean restart file 
    340340 
    341341#if defined key_agrif 
     
    361361                               CALL dia_wri_state( 'output.abort', kstp ) 
    362362      ENDIF 
     363      IF( ln_harm_ana_store   )   CALL harm_ana( kstp )        ! Harmonic analysis of tides  
    363364      IF( kstp == nit000   )   THEN 
    364365                 CALL iom_close( numror )     ! close input  ocean restart file 
     
    366367         IF( lwm.AND.numoni /= -1 ) CALL FLUSH    ( numoni )     ! flush output namelist ice 
    367368      ENDIF 
     369      IF( lrst_oce         )   CALL rst_write( kstp )       ! write output ocean restart file 
    368370 
    369371      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
  • branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/step_oce.F90

    r8058 r8059  
    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 
    109113 
  • branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/trc_oce.F90

    r8058 r8059  
    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 
Note: See TracChangeset for help on using the changeset viewer.