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

Changeset 15422


Ignore:
Timestamp:
2021-10-21T11:19:25+02:00 (2 years ago)
Author:
jcastill
Message:

Changes tested so that they can merged with the CO9 Met Office branch - jpmax_harmo should be 34 with FES14 tides, but the last components are not used anyway

Location:
NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE
Files:
2 added
15 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/BDY/bdyini.F90

    r14075 r15422  
    167167      ! Check and write out namelist parameters 
    168168      ! ----------------------------------------- 
    169       IF( jperio /= 0 )   CALL ctl_stop( 'bdy_segs: Cyclic or symmetric,',   & 
    170          &                               ' and general open boundary condition are not compatible' ) 
     169!      IF( jperio /= 0 )   CALL ctl_stop( 'bdy_segs: Cyclic or symmetric,',   & 
     170!         &                               ' and general open boundary condition are not compatible' ) 
    171171 
    172172      IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy 
  • NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DIA/diaharm.F90

    r14075 r15422  
    55   !!====================================================================== 
    66   !! History :  3.1  !  2007  (O. Le Galloudec, J. Chanut)  Original code 
     7   !!   
     8   !!   NB: 2017-12 : add 3D harmonic analysis of velocities   
     9   !!                 integration of Maria Luneva's development   
     10   !!   'key_3Ddiaharm 
    711   !!---------------------------------------------------------------------- 
    812   USE oce             ! ocean dynamics and tracers variables 
     
    1317   USE sbctide         ! Tidal forcing or not 
    1418   ! 
     19# if defined key_3Ddiaharm   
     20   USE zdf_oce   
     21#endif   
     22   ! 
    1523   USE in_out_manager  ! I/O units 
    1624   USE iom             ! I/0 library 
     
    3341   INTEGER         ::   nb_ana        ! Number of harmonics to analyse 
    3442 
    35    INTEGER , ALLOCATABLE, DIMENSION(:)       ::   name 
    36    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ana_temp 
    37    REAL(wp), ALLOCATABLE, DIMENSION(:)       ::   ana_freq, ut, vt, ft 
    38    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   out_eta, out_u, out_v 
     43   INTEGER , ALLOCATABLE, DIMENSION(:)           ::   name  
     44   REAL(wp), ALLOCATABLE, DIMENSION(:)           ::   ana_freq, ut, vt, ft   
     45# if defined key_3Ddiaharm   
     46   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:,:)   ::   ana_temp   
     47   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:)     ::   out_eta, out_u, out_v, out_w, out_dzi   
     48# else   
     49   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:)     ::   ana_temp   
     50   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)       ::   out_eta, out_u, out_v   
     51# endif 
    3952 
    4053   INTEGER  ::   ninco, nsparse 
     
    7689         WRITE(numout,*) 
    7790         WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization' 
     91# if defined key_3Ddiaharm   
     92         WRITE(numout,*) '  - 3D harmonic analysis of currents activated (key_3Ddiaharm)'   
     93#endif 
    7894         WRITE(numout,*) '~~~~~~~~~~~~~ ' 
    7995      ENDIF 
     
    155171         ! Initialize temporary arrays: 
    156172         ! ---------------------------- 
    157          ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) ) 
    158          ana_temp(:,:,:,:) = 0._wp 
     173# if defined key_3Ddiaharm   
     174         ALLOCATE( ana_temp( jpi, jpj, 2*nb_ana, 5, jpk ) )   
     175         ana_temp(:,:,:,:,:) = 0._wp   
     176# else   
     177         ALLOCATE( ana_temp( jpi, jpj, 2*nb_ana, 3      ) )   
     178         ana_temp(:,:,:,:  ) = 0._wp   
     179#endif 
    159180 
    160181      ENDIF 
     
    175196      ! 
    176197      INTEGER  ::   ji, jj, jh, jc, nhc 
     198# if defined key_3Ddiaharm   
     199      INTEGER  :: jk   
     200# endif 
    177201      REAL(wp) ::   ztime, ztemp 
    178202      !!-------------------------------------------------------------------- 
     
    190214                  &    +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh))) 
    191215                  ! 
     216               ! ssh, ub, vb are stored at the last level of 5d array 
    192217               DO jj = 2, jpjm1 
    193218                  DO ji = 2, jpim1 
    194                      ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp * sshn(ji,jj) * ssmask (ji,jj) ! elevation       
    195                      ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp * un_b(ji,jj) * ssumask(ji,jj) ! u-vel 
    196                      ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp * vn_b(ji,jj) * ssvmask(ji,jj) ! v-vel 
     219                     ! Elevation and currents   
     220# if defined key_3Ddiaharm   
     221                     ana_temp(ji,jj,nhc,1,jpk) = ana_temp(ji,jj,nhc,1,jpk) + ztemp*sshn(ji,jj)*ssmask (ji,jj)           
     222                     ana_temp(ji,jj,nhc,2,jpk) = ana_temp(ji,jj,nhc,2,jpk) + ztemp*un_b(ji,jj)*ssumask(ji,jj)   
     223                     ana_temp(ji,jj,nhc,3,jpk) = ana_temp(ji,jj,nhc,3,jpk) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj)   
     224   
     225                     ana_temp(ji,jj,nhc,5,jpk) = ana_temp(ji,jj,nhc,5,jpk)                               &   
     226                   &                              + ztemp*bfrva(ji,jj)*vn(ji,jj,mbkv(ji,jj))*ssvmask(ji,jj)   
     227                     ana_temp(ji,jj,nhc,4,jpk) = ana_temp(ji,jj,nhc,4,jpk)                               &    
     228                   &                              + ztemp*bfrua(ji,jj)*un(ji,jj,mbku(ji,jj))*ssumask(ji,jj)   
     229# else   
     230                      ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*ssmask (ji,jj)           
     231                      ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj)   
     232                      ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj)   
     233# endif 
    197234                  END DO 
    198235               END DO 
     236               !  
     237# if defined key_3Ddiaharm   
     238! 3d velocity and density:   
     239             DO jk=1,jpk-1   
     240               DO jj = 1,jpj   
     241                  DO ji = 1,jpi   
     242                     ! density and velocity   
     243                     ana_temp(ji,jj,nhc,1,jk) = ana_temp(ji,jj,nhc,1,jk) + ztemp*rhd(ji,jj,jk)   
     244                     ana_temp(ji,jj,nhc,2,jk) = ana_temp(ji,jj,nhc,2,jk) + ztemp*(un(ji,jj,jk)-un_b(ji,jj)) &   
     245                &                                          *umask(ji,jj,jk)   
     246                     ana_temp(ji,jj,nhc,3,jk) = ana_temp(ji,jj,nhc,3,jk) + ztemp*(vn(ji,jj,jk)-vn_b(ji,jj)) &   
     247                &                                          *vmask(ji,jj,jk)    
     248                     ana_temp(ji,jj,nhc,4,jk) = ana_temp(ji,jj,nhc,4,jk) + ztemp*wn(ji,jj,jk)   
     249    
     250                     ana_temp(ji,jj,nhc,5,jk) = ana_temp(ji,jj,nhc,5,jk) - 0.5*grav*ztemp*(rhd(ji,jj,jk)+rhd(ji,jj,jk+1))/max(rn2(ji,jj,jk),1.e-8_wp)   
     251                  END DO   
     252               END DO   
     253             ENDDO   
     254# endif 
    199255            END DO 
    200256         END DO 
     
    218274      !!-------------------------------------------------------------------- 
    219275      INTEGER  ::   ji, jj, jh, jc, jn, nhan 
     276# if defined key_3Ddiaharm   
     277      INTEGER :: jk   
     278# endif  
    220279      INTEGER  ::   ksp, kun, keq 
    221280      REAL(wp) ::   ztime, ztime_ini, ztime_end, z1_han 
     
    226285      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    227286       
    228       ALLOCATE( out_eta(jpi,jpj,2*nb_ana), out_u(jpi,jpj,2*nb_ana), out_v(jpi,jpj,2*nb_ana) ) 
    229  
    230287      ztime_ini = nit000_han*rdt                 ! Initial time in seconds at the beginning of analysis 
    231288      ztime_end = nitend_han*rdt                 ! Final time in seconds at the end of analysis 
     
    233290      z1_han = 1._wp / REAL(nhan-1)  
    234291       
     292# if defined key_3Ddiaharm   
     293      ALLOCATE( out_eta(jpi,jpj,jpk,2*nb_ana),   &   
     294         &      out_u  (jpi,jpj,jpk,2*nb_ana),   &   
     295         &      out_v  (jpi,jpj,jpk,2*nb_ana),   &   
     296         &      out_w  (jpi,jpj,jpk,2*nb_ana),   &   
     297         &      out_dzi(jpi,jpj,jpk,2*nb_ana) )   
     298# else   
     299      ALLOCATE( out_eta(jpi,jpj,2*nb_ana),   &   
     300         &      out_u  (jpi,jpj,2*nb_ana),   &   
     301         &      out_v  (jpi,jpj,2*nb_ana)  )   
     302# endif   
     303   
     304      IF(lwp) WRITE(numout,*) 'ANA F OLD', ft    
     305      IF(lwp) WRITE(numout,*) 'ANA U OLD', ut   
     306      IF(lwp) WRITE(numout,*) 'ANA V OLD', vt 
     307 
    235308      ninco = 2*nb_ana 
    236309 
     
    260333      CALL SUR_DETERMINE_INIT 
    261334 
    262       ! Elevation: 
     335      ! Density and Elevation:   
     336# if defined key_3Ddiaharm   
     337    DO jk=1,jpk   
     338# endif 
    263339      DO jj = 2, jpjm1 
    264340         DO ji = 2, jpim1 
    265341 
    266342            ! Fill input array 
     343# if defined key_3Ddiaharm   
     344            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,1,jk)   
     345# else 
    267346            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,1) 
     347# endif 
    268348            CALL SUR_DETERMINE 
    269349             
    270350            ! Fill output array 
    271351            DO jh = 1, nb_ana 
    272                out_eta(ji,jj,jh       ) =  ztmp7((jh-1)*2+1) * ssmask(ji,jj) 
    273                out_eta(ji,jj,jh+nb_ana) = -ztmp7((jh-1)*2+2) * ssmask(ji,jj) 
     352# if defined key_3Ddiaharm   
     353               out_eta(ji,jj,jk,jh       ) =  ztmp7((jh-1)*2+1) * ssmask(ji,jj) 
     354               out_eta(ji,jj,jk,jh+nb_ana) = -ztmp7((jh-1)*2+2) * ssmask(ji,jj) 
     355# else   
     356               out_eta(ji,jj,   jh       ) =  ztmp7((jh-1)*2+1) * ssmask(ji,jj) 
     357               out_eta(ji,jj,   jh+nb_ana) = -ztmp7((jh-1)*2+2) * ssmask(ji,jj) 
     358# endif 
    274359            END DO 
    275360         END DO 
     
    281366 
    282367            ! Fill input array 
     368# if defined key_3Ddiaharm   
     369            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,2,jk) 
     370# else  
    283371            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,2) 
     372# endif 
    284373            CALL SUR_DETERMINE 
    285374 
    286375            ! Fill output array 
    287376            DO jh = 1, nb_ana 
    288                out_u(ji,jj,       jh) =  ztmp7((jh-1)*2+1) * ssumask(ji,jj) 
    289                out_u(ji,jj,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssumask(ji,jj) 
     377# if defined key_3Ddiaharm   
     378               out_u(ji,jj,jk,       jh) =  ztmp7((jh-1)*2+1) * ssumask(ji,jj) 
     379               out_u(ji,jj,jk,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssumask(ji,jj) 
     380# else   
     381               out_u(ji,jj,          jh) =  ztmp7((jh-1)*2+1) * ssumask(ji,jj) 
     382               out_u(ji,jj,   nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssumask(ji,jj) 
     383# endif 
     384 
    290385            END DO 
    291386 
     
    298393 
    299394            ! Fill input array 
     395# if defined key_3Ddiaharm   
     396            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,3,jk) 
     397# else 
    300398            ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,3) 
     399# endif 
    301400            CALL SUR_DETERMINE 
    302401 
    303402            ! Fill output array 
    304403            DO jh = 1, nb_ana 
    305                out_v(ji,jj,       jh) =  ztmp7((jh-1)*2+1) * ssvmask(ji,jj) 
    306                out_v(ji,jj,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssvmask(ji,jj) 
     404# if defined key_3Ddiaharm   
     405               out_v(ji,jj,jk,       jh) =  ztmp7((jh-1)*2+1) * ssvmask(ji,jj) 
     406               out_v(ji,jj,jk,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssvmask(ji,jj) 
     407# else   
     408               out_v(ji,jj,          jh) =  ztmp7((jh-1)*2+1) * ssvmask(ji,jj) 
     409               out_v(ji,jj,   nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssvmask(ji,jj) 
     410# endif 
    307411            END DO 
    308412 
    309413         END DO 
    310414      END DO 
     415 
     416# if defined key_3Ddiaharm   
     417      ! w- velocity   
     418      DO jj = 1, jpj   
     419         DO ji = 1, jpi   
     420            ! Fill input array   
     421            kun=0   
     422            DO jh = 1,nb_ana   
     423               DO jc = 1,2   
     424                  kun = kun + 1   
     425                  ztmp4(kun)=ana_temp(ji,jj,kun,4,jk)   
     426               END DO   
     427            END DO   
     428   
     429            CALL SUR_DETERMINE(jj+1)   
     430   
     431            ! Fill output array   
     432            DO jh = 1, nb_ana   
     433               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)   
     434               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)   
     435            END DO   
     436   
     437         END DO   
     438      END DO   
     439   
     440      DO jj = 1, jpj   
     441         DO ji = 1, jpi   
     442            DO jh = 1, nb_ana   
     443               X1=ana_amp(ji,jj,jh,1)   
     444               X2=-ana_amp(ji,jj,jh,2)   
     445               out_w(ji,jj,jk,       jh)=X1 * tmask_i(ji,jj)   
     446               out_w(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj)   
     447            END DO   
     448         END DO   
     449      END DO   
     450   
     451       ! dzi- isopycnal displacements   
     452      DO jj = 1, jpj   
     453         DO ji = 1, jpi   
     454            ! Fill input array   
     455            kun=0   
     456            DO jh = 1,nb_ana   
     457               DO jc = 1,2   
     458                  kun = kun + 1   
     459                  ztmp4(kun)=ana_temp(ji,jj,kun,5,jk)   
     460               END DO   
     461            END DO   
     462   
     463            CALL SUR_DETERMINE(jj+1)   
     464   
     465            ! Fill output array   
     466            DO jh = 1, nb_ana   
     467               ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)   
     468               ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)   
     469            END DO   
     470   
     471         END DO   
     472      END DO   
     473   
     474      DO jj = 1, jpj   
     475         DO ji = 1, jpi   
     476            DO jh = 1, nb_ana   
     477               X1=ana_amp(ji,jj,jh,1)   
     478               X2=-ana_amp(ji,jj,jh,2)   
     479               out_dzi(ji,jj,jk,       jh)=X1 * tmask_i(ji,jj)   
     480               out_dzi(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj)   
     481            END DO   
     482         END DO   
     483      END DO   
     484   
     485   ENDDO ! jk   
     486# endif 
    311487      ! 
    312488      ! clem: we could avoid this call if all the loops were from 1:jpi and 1:jpj 
     
    328504      !!-------------------------------------------------------------------- 
    329505      INTEGER  ::   jh 
    330       !!---------------------------------------------------------------------- 
     506  
     507# if defined key_3Ddiaharm   
     508      CHARACTER(LEN=lc) :: cdfile_name_W         ! name of the file created (W-points)   
     509      INTEGER  :: jk   
     510      REAL(WP), ALLOCATABLE, DIMENSION (:,:,:) :: z3real, z3im    
     511      REAL(WP), ALLOCATABLE, DIMENSION (:,:)   :: z2real, z2im         
     512# endif   
     513!!----------------------------------------------------------------------   
     514   
     515#if defined key_dimgout   
     516      cdfile_name_T = TRIM(cexper)//'_Tidal_harmonics_gridT.dimgproc'   
     517      cdfile_name_U = TRIM(cexper)//'_Tidal_harmonics_gridU.dimgproc'   
     518      cdfile_name_V = TRIM(cexper)//'_Tidal_harmonics_gridV.dimgproc'   
     519#   if defined key_3Ddiaharm   
     520      cdfile_name_W = TRIM(cexper)//'_Tidal_harmonics_gridW.dimgproc'   
     521#   endif   
     522#endif 
    331523 
    332524      IF(lwp) WRITE(numout,*) '  ' 
    333525      IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results' 
     526#if defined key_dimgout   
     527      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~  Output files: ', TRIM(cdfile_name_T)   
     528      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_U)   
     529      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_V)   
     530#   if defined key_3Ddiaharm   
     531      IF(lwp) WRITE(numout,*) '                             ', TRIM(cdfile_name_W)   
     532#   endif   
     533#endif 
    334534      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    335535 
    336       ! A) Elevation 
     536# if defined key_3Ddiaharm   
     537      ALLOCATE(z3real(jpi,jpj,jpk),z3im(jpi,jpj,jpk),z2real(jpi,jpj),z2im(jpi,jpj))   
     538# endif   
     539   
     540      ! A) density and elevation 
    337541      !///////////// 
     542#if defined key_dimgout   
     543      cltext='density amplitude and phase; elevation is level=jpk '   
     544      CALL dia_wri_dimg(TRIM(cdfile_name_T), TRIM(cltext), out_eta, 2*nb_ana, '2')   
     545#else   
     546#   if defined key_3Ddiaharm   
     547      z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp   
     548#   endif 
    338549      DO jh = 1, nb_ana 
     550#   if defined key_3Ddiaharm   
     551        DO jk=1,jpkm1   
     552          z3real(:,:,jk)=out_eta(:,:,jk,jh)   
     553          z3im  (:,:,jk)=out_eta(:,:,jk,jh+nb_ana)   
     554        ENDDO   
     555      z2real(:,:)=out_eta(:,:,jpk,jh); z2im(:,:)=out_eta(:,:,jpk,jh+nb_ana)   
     556      CALL iom_put( TRIM(tname(jh))//'x_ro', z3real(:,:,:) )   
     557      CALL iom_put( TRIM(tname(jh))//'y_ro', z3im  (:,:,:) )   
     558      CALL iom_put( TRIM(tname(jh))//'x'   , z2real(:,:  ) )   
     559      CALL iom_put( TRIM(tname(jh))//'y'   , z2im  (:,:  ) )   
     560#   else    
     561      WRITE(numout,*) "OUTPUT ORI: ", TRIM(tname(jh))//'x', ' & ', TRIM(tname(jh))//'y', MAXVAL(out_eta(:,:,jh)) 
    339562      CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) ) 
    340563      CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,jh+nb_ana) ) 
    341       END DO 
    342  
    343       ! B) ubar 
     564#   endif   
     565      END DO   
     566#endif   
     567   
     568      ! B) u 
    344569      !///////// 
     570#if defined key_dimgout   
     571      cltext='3d u amplitude and phase; ubar is the last level'   
     572      CALL dia_wri_dimg(TRIM(cdfile_name_U), TRIM(cltext), out_u, 2*nb_ana, '2')   
     573#else   
     574#   if defined key_3Ddiaharm   
     575      z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp   
     576#   endif 
    345577      DO jh = 1, nb_ana 
    346       CALL iom_put( TRIM(tname(jh))//'x_u', out_u(:,:,jh) ) 
    347       CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,jh+nb_ana) ) 
    348       END DO 
    349  
    350       ! C) vbar 
     578#   if defined key_3Ddiaharm   
     579        DO jk=1,jpkm1   
     580          z3real(:,:,jk)=out_u(:,:,jk,jh)   
     581          z3im  (:,:,jk)=out_u(:,:,jk,jh+nb_ana)   
     582        ENDDO   
     583        z2real(:,:)=out_u(:,:,jpk,jh); z2im(:,:)=out_u(:,:,jpk,jh+nb_ana)   
     584        CALL iom_put( TRIM(tname(jh))//'x_u3d', z3real(:,:,:) )   
     585        CALL iom_put( TRIM(tname(jh))//'y_u3d', z3im (:,:,:)  )   
     586        CALL iom_put( TRIM(tname(jh))//'x_u2d', z2real(:,:) )   
     587        CALL iom_put( TRIM(tname(jh))//'y_u2d', z2im (:,:)  )   
     588        z2real(:,:)=out_w(:,:,jpk,jh); z2im(:,:)=out_w(:,:,jpk,jh+nb_ana)   
     589        CALL iom_put( TRIM(tname(jh))//'x_tabx', z2real(:,:) )   
     590        CALL iom_put( TRIM(tname(jh))//'y_tabx', z2im (:,:)  )   
     591#   else   
     592        CALL iom_put( TRIM(tname(jh))//'x_u2d', out_u(:,:,jh) )   
     593        CALL iom_put( TRIM(tname(jh))//'y_u2d', out_u(:,:,nb_ana+jh) )   
     594#   endif   
     595      END DO   
     596#endif   
     597   
     598      ! C) v 
    351599      !///////// 
     600#if defined key_dimgout   
     601      cltext='3d v amplitude and phase; vbar is the last level'   
     602      CALL dia_wri_dimg(TRIM(cdfile_name_V), TRIM(cltext), out_v, 2*nb_ana, '2')   
     603#else   
     604#   if defined key_3Ddiaharm   
     605      z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp   
     606#   endif 
    352607      DO jh = 1, nb_ana 
    353          CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh       ) ) 
    354          CALL iom_put( TRIM(tname(jh))//'y_v', out_v(:,:,jh+nb_ana) ) 
    355       END DO 
    356       ! 
     608#   if defined key_3Ddiaharm   
     609        DO jk=1,jpkm1   
     610          z3real(:,:,jk)=out_v(:,:,jk,jh)   
     611          z3im  (:,:,jk)=out_v(:,:,jk,jh+nb_ana)   
     612        ENDDO   
     613        z2real(:,:)=out_v(:,:,jpk,jh); z2im(:,:)=out_v(:,:,jpk,jh+nb_ana)   
     614        CALL iom_put( TRIM(tname(jh))//'x_v3d', z3real(:,:,:) )   
     615        CALL iom_put( TRIM(tname(jh))//'y_v3d', z3im (:,:,:)  )   
     616        CALL iom_put( TRIM(tname(jh))//'x_v2d'  , z2real(:,:) )   
     617        CALL iom_put( TRIM(tname(jh))//'y_v2d'  , z2im (:,:)  )   
     618        z2real(:,:)=out_dzi(:,:,jpk,jh); z2im(:,:)=out_dzi(:,:,jpk,jh+nb_ana)   
     619        CALL iom_put( TRIM(tname(jh))//'x_taby', z2real(:,:) )   
     620        CALL iom_put( TRIM(tname(jh))//'y_taby', z2im (:,:)  )   
     621#   else   
     622         CALL iom_put( TRIM(tname(jh))//'x_v2d', out_v(:,:,jh ) )   
     623         CALL iom_put( TRIM(tname(jh))//'y_v2d', out_v(:,:,jh+nb_ana) )   
     624#   endif   
     625       END DO   
     626   
     627#endif   
     628      ! D) w   
     629# if defined key_3Ddiaharm   
     630#   if defined key_dimgout   
     631      cltext='3d w amplitude and phase; vort_baro is the last level'   
     632      CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2')   
     633#   else   
     634      DO jh = 1, nb_ana   
     635        DO jk=1,jpkm1   
     636         z3real(:,:,jk)=out_w(:,:,jk,jh)   
     637         z3im(:,:,jk)=out_w(:,:,jk,jh+nb_ana)   
     638        ENDDO   
     639        CALL iom_put( TRIM(tname(jh))//'x_w3d', z3real(:,:,:) )   
     640        CALL iom_put( TRIM(tname(jh))//'y_w3d', z3im(:,:,:) )   
     641      END DO   
     642#   endif   
     643   
     644!       E) dzi + tau_bot   
     645#   if defined key_dimgout   
     646      cltext='dzi=g*ro/N2 amplitude and phase'   
     647      CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2')   
     648#   else   
     649      DO jh = 1, nb_ana   
     650        DO jk=1,jpkm1   
     651         z3real(:,:,jk)=out_dzi(:,:,jk,jh)   
     652         z3im(:,:,jk)=out_dzi(:,:,jk,jh+nb_ana)   
     653        ENDDO   
     654        CALL iom_put( TRIM(tname(jh))//'x_dzi', z3real(:,:,:) )   
     655        CALL iom_put( TRIM(tname(jh))//'y_dzi', z3im(:,:,:) )   
     656      END DO   
     657#   endif   
     658# endif    
     659   
     660      !   
     661# if defined key_3Ddiaharm   
     662   DEALLOCATE(z3real, z3im, z2real,z2im)   
     663# endif 
     664 
    357665   END SUBROUTINE dia_wri_harm 
    358666 
  • NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DOM/dtatsd.F90

    r14075 r15422  
    2121   ! 
    2222   USE in_out_manager  ! I/O manager 
     23   USE iom             ! IOM library 
    2324   USE lib_mpp         ! MPP library 
    2425 
     
    3132   !                                  !!* namtsd  namelist : Temperature & Salinity Data * 
    3233   LOGICAL , PUBLIC ::   ln_tsd_init   !: T & S data flag 
     34   LOGICAL , PUBLIC ::   ln_tsd_interp !: vertical interpolation flag 
    3335   LOGICAL , PUBLIC ::   ln_tsd_dmp    !: internal damping toward input data flag 
    3436 
    3537   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_tsd   ! structure of input SST (file informations, fields read) 
     38   INTEGER                              ::   jpk_init, inum_dta   
     39   INTEGER                              ::   id, linum   ! local integers   
     40   INTEGER                              ::   zdim(4) 
    3641 
    3742   !!---------------------------------------------------------------------- 
     
    5358      LOGICAL, INTENT(in), OPTIONAL ::   ld_tradmp   ! force the initialization when tradp is used 
    5459      ! 
    55       INTEGER ::   ios, ierr0, ierr1, ierr2, ierr3   ! local integers 
     60      INTEGER ::   ios, ierr0, ierr1, ierr2, ierr3, ierr4, ierr5   ! local integers 
    5661      !! 
    5762      CHARACTER(len=100)            ::   cn_dir          ! Root directory for location of ssr files 
    58       TYPE(FLD_N), DIMENSION( jpts) ::   slf_i           ! array of namelist informations on the fields to read 
    59       TYPE(FLD_N)                   ::   sn_tem, sn_sal 
     63      TYPE(FLD_N), DIMENSION(jpts+2)::   slf_i           ! array of namelist informations on the fields to read   
     64      TYPE(FLD_N)                   ::   sn_tem, sn_sal, sn_dep, sn_msk  
    6065      !! 
    6166      NAMELIST/namtsd/   ln_tsd_init, ln_tsd_dmp, cn_dir, sn_tem, sn_sal 
     
    6368      ! 
    6469      !  Initialisation 
    65       ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0 
     70      ierr0 = 0  ;  ierr1 = 0  ;  ierr2 = 0  ;  ierr3 = 0  ; ierr4 = 0  ;  ierr5 = 0 
    6671      ! 
    6772      REWIND( numnam_ref )              ! Namelist namtsd in reference namelist :  
     
    8085         WRITE(numout,*) '~~~~~~~~~~~~ ' 
    8186         WRITE(numout,*) '   Namelist namtsd' 
    82          WRITE(numout,*) '      Initialisation of ocean T & S with T &S input data   ln_tsd_init = ', ln_tsd_init 
    83          WRITE(numout,*) '      damping of ocean T & S toward T &S input data        ln_tsd_dmp  = ', ln_tsd_dmp 
     87         WRITE(numout,*) '      Initialisation of ocean T & S with T & S input data  ln_tsd_init   = ', ln_tsd_init  
     88         WRITE(numout,*) '      Interpolation of initial conditions in the vertical  ln_tsd_interp = ', ln_tsd_interp  
     89         WRITE(numout,*) '      damping of ocean T & S toward T & S input data       ln_tsd_dmp    = ', ln_tsd_dmp 
    8490         WRITE(numout,*) 
    8591         IF( .NOT.ln_tsd_init .AND. .NOT.ln_tsd_dmp ) THEN 
     
    94100         ln_tsd_init = .FALSE. 
    95101      ENDIF 
     102      IF( ln_tsd_interp .AND. ln_tsd_dmp ) THEN   
     103            CALL ctl_stop( 'dta_tsd_init: Tracer damping and vertical interpolation not yet configured' )   ;   RETURN   
     104      ENDIF   
     105      IF( ln_tsd_interp .AND. LEN(TRIM(sn_msk%wname)) > 0 ) THEN   
     106            CALL ctl_stop( 'dta_tsd_init: Using vertical interpolation and weights files not recommended' )   ;   RETURN   
     107      ENDIF 
    96108      ! 
    97109      !                             ! allocate the arrays (if necessary) 
    98110      IF( ln_tsd_init .OR. ln_tsd_dmp ) THEN 
    99111         ! 
    100          ALLOCATE( sf_tsd(jpts), STAT=ierr0 ) 
     112         IF( ln_tsd_interp ) THEN   
     113           ALLOCATE( sf_tsd(jpts+2), STAT=ierr0 ) ! to carry the addtional depth information   
     114         ELSE   
     115           ALLOCATE( sf_tsd(jpts  ), STAT=ierr0 )    
     116         ENDIF 
    101117         IF( ierr0 > 0 ) THEN 
    102118            CALL ctl_stop( 'dta_tsd_init: unable to allocate sf_tsd structure' )   ;   RETURN 
    103119         ENDIF 
    104120         ! 
    105                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk)   , STAT=ierr0 ) 
    106          IF( sn_tem%ln_tint )   ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 ) 
    107                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk)   , STAT=ierr2 ) 
    108          IF( sn_sal%ln_tint )   ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 ) 
    109          ! 
    110          IF( ierr0 + ierr1 + ierr2 + ierr3 > 0 ) THEN 
     121         IF( ln_tsd_interp ) THEN   
     122            CALL iom_open ( trim(cn_dir) // trim(sn_dep%clname), inum_dta )    
     123            id = iom_varid( inum_dta, sn_dep%clvar, zdim )   
     124            jpk_init = zdim(3)   
     125            IF(lwp) WRITE(numout,*) 'Dimension of vertical coordinate in ICs: ', jpk_init   
     126            CALL iom_close( inum_dta )   ! Close the input file   
     127            !   
     128                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr0 )   
     129            IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr1 )   
     130                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr2 )   
     131            IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr3 )     
     132                                 ALLOCATE( sf_tsd(jp_dep)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr4 )   
     133                                 ALLOCATE( sf_tsd(jp_msk)%fnow(jpi,jpj,jpk_init  ) , STAT=ierr5 )   
     134         ELSE   
     135                                 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk)   , STAT=ierr0 )   
     136            IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 )   
     137                                 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk)   , STAT=ierr2 )   
     138            IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 )     
     139         ENDIF ! ln_tsd_interp   
     140         !   
     141         IF( ierr0 + ierr1 + ierr2 + ierr3 + ierr4 + ierr5 > 0 ) THEN  
    111142            CALL ctl_stop( 'dta_tsd : unable to allocate T & S data arrays' )   ;   RETURN 
    112143         ENDIF 
    113144         !                         ! fill sf_tsd with sn_tem & sn_sal and control print 
    114145         slf_i(jp_tem) = sn_tem   ;   slf_i(jp_sal) = sn_sal 
     146         IF( ln_tsd_interp ) slf_i(jp_dep) = sn_dep   ;   slf_i(jp_msk) = sn_msk 
    115147         CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd', no_print ) 
    116148         ! 
     
    138170      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(  out) ::   ptsd   ! T & S data 
    139171      ! 
    140       INTEGER ::   ji, jj, jk, jl, jkk   ! dummy loop indicies 
    141       INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1   ! local integers 
     172      INTEGER ::   ji, jj, jk, jl, jk_init   ! dummy loop indices   
     173      INTEGER ::   ik, il0, il1, ii0, ii1, ij0, ij1        ! local integers  
    142174      REAL(wp)::   zl, zi                             ! local scalars 
    143       REAL(wp), DIMENSION(jpk) ::  ztp, zsp   ! 1D workspace 
    144175      !!---------------------------------------------------------------------- 
    145176      ! 
     
    176207!!gm end 
    177208      ! 
    178       ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:)    ! NO mask 
    179       ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:)  
    180       ! 
    181       IF( ln_sco ) THEN                   !==   s- or mixed s-zps-coordinate   ==! 
    182          ! 
    183          IF( kt == nit000 .AND. lwp )THEN 
    184             WRITE(numout,*) 
    185             WRITE(numout,*) 'dta_tsd: interpolates T & S data onto the s- or mixed s-z-coordinate mesh' 
    186          ENDIF 
    187          ! 
    188          DO jj = 1, jpj                         ! vertical interpolation of T & S 
    189             DO ji = 1, jpi 
    190                DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points 
     209      IF( kt == nit000 .AND. lwp )THEN   
     210         WRITE(numout,*)   
     211         WRITE(numout,*) 'dta_tsd: interpolates T & S data onto current mesh'   
     212      ENDIF   
     213      !   
     214      IF( ln_tsd_interp ) THEN                 ! probably should use pointers in the following to make more readable   
     215      !   
     216         DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points   
     217            DO jj= 1, jpj   
     218               DO ji= 1, jpi 
    191219                  zl = gdept_0(ji,jj,jk) 
    192                   IF(     zl < gdept_1d(1  ) ) THEN          ! above the first level of data 
    193                      ztp(jk) =  ptsd(ji,jj,1    ,jp_tem) 
    194                      zsp(jk) =  ptsd(ji,jj,1    ,jp_sal) 
    195                   ELSEIF( zl > gdept_1d(jpk) ) THEN          ! below the last level of data 
    196                      ztp(jk) =  ptsd(ji,jj,jpkm1,jp_tem) 
    197                      zsp(jk) =  ptsd(ji,jj,jpkm1,jp_sal) 
    198                   ELSE                                      ! inbetween : vertical interpolation between jkk & jkk+1 
    199                      DO jkk = 1, jpkm1                                  ! when  gdept(jkk) < zl < gdept(jkk+1) 
    200                         IF( (zl-gdept_1d(jkk)) * (zl-gdept_1d(jkk+1)) <= 0._wp ) THEN 
    201                            zi = ( zl - gdept_1d(jkk) ) / (gdept_1d(jkk+1)-gdept_1d(jkk)) 
    202                            ztp(jk) = ptsd(ji,jj,jkk,jp_tem) + ( ptsd(ji,jj,jkk+1,jp_tem) - ptsd(ji,jj,jkk,jp_tem) ) * zi  
    203                            zsp(jk) = ptsd(ji,jj,jkk,jp_sal) + ( ptsd(ji,jj,jkk+1,jp_sal) - ptsd(ji,jj,jkk,jp_sal) ) * zi 
     220                  IF( zl < sf_tsd(jp_dep)%fnow(ji,jj,1) ) THEN                     ! above the first level of data   
     221                     ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,1)    
     222                     ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,1)   
     223                  ELSEIF( zl > sf_tsd(jp_dep)%fnow(ji,jj,jpk_init) ) THEN          ! below the last level of data   
     224                     ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jpk_init)   
     225                     ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jpk_init)   
     226                  ELSE                                                             ! inbetween : vertical interpolation between jk_init & jk_init+1   
     227                     DO jk_init = 1, jpk_init-1                                    ! when  gdept(jk_init) < zl < gdept(jk_init+1)   
     228                        IF( sf_tsd(jp_msk)%fnow(ji,jj,jk_init+1) == 0 ) THEN       ! if there is no data fill down   
     229                           sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init)   
     230                           sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init)   
     231                        ENDIF   
     232                        IF((zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init)) * (zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)) <= 0._wp ) THEN   
     233                           zi = ( zl - sf_tsd(jp_dep)%fnow(ji,jj,jk_init) ) / &   
     234                        &       (sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_dep)%fnow(ji,jj,jk_init))   
     235                           ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init) + &   
     236                        &                         (sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_tem)%fnow(ji,jj,jk_init)) * zi   
     237                           ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init) + &   
     238                        &                         (sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_sal)%fnow(ji,jj,jk_init)) * zi 
    204239                        ENDIF 
    205240                     END DO 
    206241                  ENDIF 
    207242               END DO 
    208                DO jk = 1, jpkm1 
    209                   ptsd(ji,jj,jk,jp_tem) = ztp(jk) * tmask(ji,jj,jk)     ! mask required for mixed zps-s-coord 
    210                   ptsd(ji,jj,jk,jp_sal) = zsp(jk) * tmask(ji,jj,jk) 
    211                END DO 
    212                ptsd(ji,jj,jpk,jp_tem) = 0._wp 
    213                ptsd(ji,jj,jpk,jp_sal) = 0._wp 
    214243            END DO 
    215244         END DO 
    216245         !  
     246         ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) *tmask(:,:,:)   
     247         ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) *tmask(:,:,:) 
    217248      ELSE                                !==   z- or zps- coordinate   ==! 
    218249         !                              
    219          ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) * tmask(:,:,:)    ! Mask 
    220          ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) * tmask(:,:,:) 
     250         ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:) * tmask(:,:,:)  ! Mask   
     251         ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:) * tmask(:,:,:) 
    221252         ! 
    222253         IF( ln_zps ) THEN                      ! zps-coordinate (partial steps) interpolation at the last ocean level 
     
    248279                                        DEALLOCATE( sf_tsd(jp_sal)%fnow )     ! S arrays in the structure 
    249280         IF( sf_tsd(jp_sal)%ln_tint )   DEALLOCATE( sf_tsd(jp_sal)%fdta ) 
     281         IF( ln_tsd_interp )            DEALLOCATE( sf_tsd(jp_dep)%fnow )     ! T arrays in the structure   
     282         IF( ln_tsd_interp )            DEALLOCATE( sf_tsd(jp_msk)%fnow )     ! T arrays in the structure 
    250283                                        DEALLOCATE( sf_tsd              )     ! the structure itself 
    251284      ENDIF 
  • NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DYN/dynnxt.F90

    r14075 r15422  
    3333   USE dynadv         ! dynamics: vector invariant versus flux form 
    3434   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme 
     35   USE dynspg 
    3536   USE domvvl         ! variable volume 
    3637   USE bdy_oce , ONLY : ln_bdy 
     
    5859   PUBLIC    dyn_nxt   ! routine called by step.F90 
    5960 
     61   !! Substitution  
     62#  include "vectopt_loop_substitute.h90" 
    6063   !!---------------------------------------------------------------------- 
    6164   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    182185            vn(:,:,jk) = va(:,:,jk) 
    183186         END DO 
     187         ! limit velocities   
     188         IF (ln_ulimit) THEN   
     189            call dyn_limit_velocity (kt)   
     190         ENDIF 
    184191         IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n 
    185192!!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a   
     
    214221               END DO 
    215222            END DO 
     223            ! limit velocities   
     224            IF (ln_ulimit) THEN   
     225               call dyn_limit_velocity (kt)   
     226            ENDIF 
    216227            !                             ! ================! 
    217228         ELSE                             ! Variable volume ! 
     
    263274                  END DO 
    264275               END DO 
     276               ! limit velocities   
     277               IF (ln_ulimit) THEN   
     278                  call dyn_limit_velocity (kt)   
     279               ENDIF 
    265280               ! 
    266281            ELSE                          ! Asselin filter applied on thickness weighted velocity 
     
    290305                  END DO 
    291306               END DO 
     307               ! limit velocities   
     308               IF (ln_ulimit) THEN   
     309                  call dyn_limit_velocity (kt)   
     310               ENDIF 
    292311               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor 
    293312               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 
     
    401420   END SUBROUTINE dyn_nxt 
    402421 
     422   SUBROUTINE dyn_limit_velocity (kt)   
     423      ! limits maximum values of un and vn by volume courant number   
     424      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index   
     425      !   
     426      INTEGER  ::   ji, jj, jk   ! dummy loop indices   
     427      REAL(wp) ::   zzu,zplim,zmlim,isp,ism,zcn,ze3e1,zzcn,zcnn,idivp,idivm   
     428    
     429      ! limit fluxes   
     430      zcn =cn_ulimit !0.9 ! maximum velocity inverse courant number   
     431      zcnn = cnn_ulimit !0.54 ! how much to reduce cn by in divergen flow   
     432    
     433      DO jk = 1, jpkm1   
     434         DO jj = 2, jpjm1   
     435            DO ji = fs_2, fs_jpim1   ! vect. opt.  
     436               ! U direction   
     437               zzu = un(ji,jj,jk)   
     438               ze3e1 = e3u_n(ji  ,jj,jk) * e2u(ji  ,jj)   
     439               ! ips is 1 if flow is positive othersize zero   
     440               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp )   
     441               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp )   
     442               ! idev is 1 if divergent flow otherwise zero   
     443               idivp = -isp * 0.5 * (sign(1.0_wp, un(ji-1,jj,jk)) - 1.0_wp )   
     444               idivm =  ism * 0.5 * (sign(1.0_wp, un(ji+1,jj,jk)) + 1.0_wp )   
     445               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp   
     446               zzcn = zzcn * zcn   
     447               zplim =  zzcn * (e3t_n(ji  ,jj,jk) * e1t(ji  ,jj) * e2t(ji  ,jj)) / &   
     448                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk)   
     449               zmlim = -zzcn * (e3t_n(ji+1,jj,jk) * e1t(ji+1,jj) * e2t(ji+1,jj)) / &   
     450                                               (2.0*rdt * ze3e1)*umask(ji,jj,jk)   
     451               ! limit currents   
     452               un(ji,jj,jk) = min ( zzu,zplim) * isp + max(zzu,zmlim) *ism   
     453    
     454               ! V  direction   
     455               zzu = vn(ji,jj,jk)   
     456               ze3e1 = e3v_n(ji  ,jj,jk) * e1v(ji  ,jj)   
     457               isp =  0.5 * (sign(1.0_wp,zzu) + 1.0_wp )   
     458               ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp )   
     459               ! idev is 1 if divergent flow otherwise zero   
     460               idivp = -isp * 0.5 * (sign(1.0_wp, vn(ji,jj-1,jk)) - 1.0_wp )   
     461               idivm =  ism * 0.5 * (sign(1.0_wp, vn(ji,jj+1,jk)) + 1.0_wp )   
     462               zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp   
     463               zzcn = zzcn * zcn   
     464               zplim =  zzcn * (e3t_n(ji,jj  ,jk) * e1t(ji,jj  ) * e2t(ji,jj  )) / &   
     465                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk)   
     466               zmlim = -zzcn * (e3t_n(ji,jj+1,jk) * e1t(ji,jj+1) * e2t(ji,jj+1)) / &   
     467                                               (2.0*rdt * ze3e1)*vmask(ji,jj,jk)   
     468               ! limit currents   
     469               vn(ji,jj,jk) = min ( zzu,zplim) * isp + max(zzu,zmlim) *ism   
     470            ENDDO   
     471         ENDDO   
     472      ENDDO   
     473       CALL lbc_lnk_multi( 'dynnxt', un(:,:,:), 'U',  -1., vn(:,:,:), 'V',  -1. )  
     474    
     475    END SUBROUTINE dyn_limit_velocity 
     476 
    403477   !!========================================================================= 
    404478END MODULE dynnxt 
  • NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DYN/dynspg.F90

    r14075 r15422  
    3939   INTEGER ::   nspg = 0   ! type of surface pressure gradient scheme defined from lk_dynspg_...  
    4040 
     41   LOGICAL, PUBLIC ::  ln_ulimit   
     42   REAL(wp), PUBLIC :: cn_ulimit,cnn_ulimit 
     43 
    4144   !                       ! Parameter to control the surface pressure gradient scheme 
    4245   INTEGER, PARAMETER ::   np_TS  = 1   ! split-explicit time stepping (Time-Splitting) 
     
    191194      NAMELIST/namdyn_spg/ ln_dynspg_exp       , ln_dynspg_ts,   & 
    192195      &                    ln_bt_fw, ln_bt_av  , ln_bt_auto  ,   & 
    193       &                    nn_baro , rn_bt_cmax, nn_bt_flt, rn_bt_alpha 
     196      &                    nn_baro , rn_bt_cmax, nn_bt_flt, rn_bt_alpha, &  
     197      &                    ln_ulimit, cn_ulimit, cnn_ulimit 
    194198      !!---------------------------------------------------------------------- 
    195199      ! 
     
    213217         WRITE(numout,*) '      Explicit free surface                  ln_dynspg_exp = ', ln_dynspg_exp 
    214218         WRITE(numout,*) '      Free surface with time splitting       ln_dynspg_ts  = ', ln_dynspg_ts 
     219         WRITE(numout,*) '     Limit velocities ln_ulimit     = ', ln_ulimit   
     220         WRITE(numout,*) '     Limit velocities      max inverse Courant number     = ', cn_ulimit   
     221         WRITE(numout,*) '     Limit velocities   multiplier for divergant flow     = ', cnn_ulimit 
    215222      ENDIF 
    216223      !                          ! Control of surface pressure gradient scheme options 
  • NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/fldread.F90

    r14075 r15422  
    901901            IF(     zdepth(jk) <= pdta_read_z(jb,1,1)      ) THEN   ! above the first level of external data 
    902902               pdta(jb,1,jk) = pdta_read(jb,1,1) 
    903             ELSEIF( zdepth(jk) >  pdta_read_z(jb,1,ipkmax) ) THEN   ! below the last level of external data /= FillValue 
    904                pdta(jb,1,jk) = pdta_read(jb,1,ipkmax) 
     903            ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkmax) ) THEN    ! below the last level of external data   
     904               pdta(jb,1,jk) =  pdta_read(jb,1,ipkmax) 
    905905            ELSE                                                    ! inbetween: vertical interpolation between jkb & jkb+1 
    906906               DO jkb = 1, ipkmax-1 
  • NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/sbctide.F90

    r14075 r15422  
    1616   USE ioipsl         ! NetCDF IPSL library 
    1717   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     18   USE bdytides 
    1819 
    1920   IMPLICIT NONE 
     
    108109 
    109110      DO jk = 1, nb_harmo 
    110          zcons = 0.7_wp * Wave(ntide(jk))%equitide * ftide(jk) 
     111         ! love number now provides in tide namelist   
     112         zcons = dn_love_number * Wave(ntide(jk))%equitide * ftide(jk) 
    111113         DO ji = 1, jpi 
    112114            DO jj = 1, jpj 
     
    119121               IF    ( Wave(ntide(jk))%nutide == 1 )  THEN  ;  zcs = zcons * SIN( 2._wp*zlat ) 
    120122               ELSEIF( Wave(ntide(jk))%nutide == 2 )  THEN  ;  zcs = zcons * COS( zlat )**2 
     123               ! Add tide potential for long period tides   
     124               ELSEIF( Wave(ntide(jk))%nutide == 0 )  THEN  ;  zcs = zcons * (0.5_wp-1.5_wp*SIN(zlat)**2._wp) 
    121125               ELSE                                         ;  zcs = 0._wp 
    122126               ENDIF 
  • NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/tide_mod.F90

    r14075 r15422  
    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 
     
    4141 
    4242   SUBROUTINE tide_init_Wave 
     43# if defined key_FES14_tides   
     44#     include "tide_FES14.h90"   
     45# else 
    4346#     include "tide.h90" 
     47# endif 
    4448   END SUBROUTINE tide_init_Wave 
    4549 
     
    331335         zf = zf * zf1 * zf1 
    332336         ! 
     337      CASE( 20 )                 !==  formule 20,  compound waves ( 78 x 78 x 78 x 78 )   
     338         zf1 = nodal_factort(78)   
     339         zf  = zf1 * zf1 * zf1 * zf1   
     340         !  
    333341      CASE( 73 )                 !==  formule 73 
    334342         zs = sin(sh_I) 
  • NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/tideini.F90

    r14075 r15422  
    3333   INTEGER , PUBLIC ::   kt_tide         !: 
    3434   REAL(wp), PUBLIC ::   rdttideramp     !: 
     35   REAL(wp), PUBLIC ::   dn_love_number  !: 
    3536   REAL(wp), PUBLIC ::   rn_scal_load    !: 
    3637   CHARACTER(lc), PUBLIC ::   cn_tide_load   !:  
     
    5455      ! 
    5556      NAMELIST/nam_tide/ln_tide, ln_tide_pot, ln_scal_load, ln_read_load, cn_tide_load, & 
    56                   &     ln_tide_ramp, rn_scal_load, rdttideramp, clname 
     57                  &     ln_tide_ramp, rn_scal_load, rdttideramp, dn_love_number, clname 
    5758      !!---------------------------------------------------------------------- 
    5859      ! 
     
    7879            WRITE(numout,*) '         Read load potential from file           ln_read_load = ', ln_read_load 
    7980            WRITE(numout,*) '         Apply ramp on tides at startup          ln_tide_ramp = ', ln_tide_ramp 
     81            WRITE(numout,*) '                                              dn_love_number  = ', dn_love_number 
    8082            WRITE(numout,*) '         Fraction of SSH used in scal. approx.   rn_scal_load = ', rn_scal_load 
    8183            WRITE(numout,*) '         Duration (days) of ramp                 rdttideramp  = ', rdttideramp 
     
    99101      END DO 
    100102      !        
     103      IF (ln_tide .and.lwp) WRITE(numout,*) 'nb_harmo     = ', nb_harmo 
     104 
    101105      ! Ensure that tidal components have been set in namelist_cfg 
    102106      IF( nb_harmo == 0 )   CALL ctl_stop( 'tide_init : No tidal components set in nam_tide' ) 
  • NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/USR/usrdef_istate.F90

    r14075 r15422  
    6262         DO jj = 1, jpj 
    6363            DO ji = 1, jpi 
    64                pts(ji,jj,jk,jp_tem) =  (  (  16. - 12. * TANH( (pdept(ji,jj,jk) - 400) / 700 ) )   & 
    65                     &           * (-TANH( (500. - pdept(ji,jj,jk)) / 150. ) + 1.) / 2.             & 
    66                     &           + ( 15. * ( 1. - TANH( (pdept(ji,jj,jk)-50.) / 1500.) )            & 
    67                     &           - 1.4 * TANH((pdept(ji,jj,jk)-100.) / 100.)                        & 
    68                     &           + 7.  * (1500. - pdept(ji,jj,jk) ) / 1500.)                        & 
    69                     &           * (-TANH( (pdept(ji,jj,jk) - 500.) / 150.) + 1.) / 2.  ) * ptmask(ji,jj,jk) 
    70  
    71                pts(ji,jj,jk,jp_sal) =  (  (  36.25 - 1.13 * TANH( (pdept(ji,jj,jk) - 305) / 460 ) )  & 
    72                     &         * (-TANH((500. - pdept(ji,jj,jk)) / 150.) + 1.) / 2                  & 
    73                     &         + ( 35.55 + 1.25 * (5000. - pdept(ji,jj,jk)) / 5000.                 & 
    74                     &         - 1.62 * TANH( (pdept(ji,jj,jk) - 60.  ) / 650. )                    & 
    75                     &         + 0.2  * TANH( (pdept(ji,jj,jk) - 35.  ) / 100. )                    & 
    76                     &         + 0.2  * TANH( (pdept(ji,jj,jk) - 1000.) / 5000.) )                  & 
    77                     &         * (-TANH( (pdept(ji,jj,jk) - 500.) / 150.) + 1.) / 2  ) * ptmask(ji,jj,jk) 
     64               pts(ji,jj,jk,jp_tem) = 20._wp * ptmask(ji,jj,jk)   
     65               pts(ji,jj,jk,jp_sal) = 36.25_wp * ptmask(ji,jj,jk) 
    7866            END DO 
    7967         END DO 
  • NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/USR/usrdef_sbc.F90

    r14075 r15422  
    11MODULE usrdef_sbc 
    22   !!====================================================================== 
    3    !!                     ***  MODULE  usrdef_sbc  *** 
    4    !! 
    5    !!                     ===  GYRE configuration  === 
     3   !!                       ***  MODULE usrdef_sbc  ***  
     4   !!   
     5   !!                  ===  AMM7_SURGE configuration  === 
    66   !! 
    77   !! User defined :   surface forcing of a user configuration 
    88   !!====================================================================== 
    99   !! History :  4.0   ! 2016-03  (S. Flavoni, G. Madec)  user defined interface 
     10   !!            4.0   ! 2017-12  (C. O'Neill)  add necessary options for surge work - either no fluxes   
     11   !!                                           (for tide-only run) or wind and pressure only 
    1012   !!---------------------------------------------------------------------- 
    1113 
    1214   !!---------------------------------------------------------------------- 
    13    !!   usrdef_sbc    : user defined surface bounday conditions in GYRE case 
     15   !!   usrdef_sbc    : user defined surface bounday conditions in LOCK_EXCHANGE case 
    1416   !!---------------------------------------------------------------------- 
    15    USE oce            ! ocean dynamics and tracers 
    16    USE dom_oce        ! ocean space and time domain 
    17    USE sbc_oce        ! Surface boundary condition: ocean fields 
    18    USE phycst         ! physical constants 
     17   USE oce             ! ocean dynamics and tracers  
     18   USE dom_oce         ! ocean space and time domain  
     19   USE sbc_oce         ! Surface boundary condition: ocean fields  
     20   USE sbc_ice         ! Surface boundary condition: ocean fields  
     21   USE fldread         ! read input fields  
     22   USE phycst          ! physical constants  
     23   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    1924   ! 
    20    USE in_out_manager ! I/O manager 
    21    USE lib_mpp        ! distribued memory computing library 
    22    USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    23    USE lib_fortran    ! 
     25   USE in_out_manager  ! I/O manager  
     26   USE iom  
     27   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)  
     28   USE lib_mpp         ! distribued memory computing library  
     29   !USE wrk_nemo       ! work arrays  
     30   USE timing         ! Timing  
     31   USE prtctl         ! Print control 
    2432 
    2533   IMPLICIT NONE 
    2634   PRIVATE 
    2735 
    28    PUBLIC   usrdef_sbc_oce       ! routine called in sbcmod module 
    29    PUBLIC   usrdef_sbc_ice_tau   ! routine called by icestp.F90 for ice dynamics 
    30    PUBLIC   usrdef_sbc_ice_flx   ! routine called by icestp.F90 for ice thermo 
     36   PUBLIC   usrdef_sbc_oce    ! routine called in sbcmod module  
     37   PUBLIC   usrdef_sbc_ice_tau  ! routine called by sbcice_lim.F90 for ice dynamics  
     38   PUBLIC   usrdef_sbc_ice_flx  ! routine called by sbcice_lim.F90 for ice thermo  
     39   !                                  !!* Namelist namsbc_usr  
     40   REAL(wp) ::   rn_vfac     ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem)  
     41   REAL(wp) ::   rn_charn_const   
     42   LOGICAL  ::   ln_use_sbc  ! Surface fluxes on or not 
    3143 
    3244   !! * Substitutions 
     
    4355      !!                    ***  ROUTINE usrdef_sbc  *** 
    4456      !!               
    45       !! ** Purpose :   provide at each time-step the GYRE surface boundary 
     57      !! ** Purpose :   provide at each time-step the surface boundary 
    4658      !!              condition, i.e. the momentum, heat and freshwater fluxes. 
    4759      !! 
    48       !! ** Method  :   analytical seasonal cycle for GYRE configuration. 
     60      !! ** Method  :   all 0 fields, for AMM7_SURGE case 
    4961      !!                CAUTION : never mask the surface stress field ! 
    5062      !! 
    51       !! ** Action  : - set the ocean surface boundary condition, i.e.    
     63      !! ** Action  : - if tide-only case - set to ZERO all the ocean surface boundary condition, i.e. 
    5264      !!                   utau, vtau, taum, wndm, qns, qsr, emp, sfx 
    53       !! 
    54       !! Reference : Hazeleger, W., and S. Drijfhout, JPO, 30, 677-695, 2000. 
     65      !!              - if tide+surge case - read in wind and air pressure      !! 
    5566      !!---------------------------------------------------------------------- 
    5667      INTEGER, INTENT(in) ::   kt   ! ocean time step 
    57       !! 
    58       INTEGER  ::   ji, jj                 ! dummy loop indices 
    59       INTEGER  ::   zyear0                 ! initial year  
    60       INTEGER  ::   zmonth0                ! initial month 
    61       INTEGER  ::   zday0                  ! initial day 
    62       INTEGER  ::   zday_year0             ! initial day since january 1st 
    63       REAL(wp) ::   ztau     , ztau_sais   ! wind intensity and of the seasonal cycle 
    64       REAL(wp) ::   ztime                  ! time in hour 
    65       REAL(wp) ::   ztimemax , ztimemin    ! 21th June, and 21th decem. if date0 = 1st january 
    66       REAL(wp) ::   ztimemax1, ztimemin1   ! 21th June, and 21th decem. if date0 = 1st january 
    67       REAL(wp) ::   ztimemax2, ztimemin2   ! 21th June, and 21th decem. if date0 = 1st january 
    68       REAL(wp) ::   ztaun                  ! intensity 
    69       REAL(wp) ::   zemp_s, zemp_n, zemp_sais, ztstar 
    70       REAL(wp) ::   zcos_sais1, zcos_sais2, ztrp, zconv, t_star 
    71       REAL(wp) ::   zsumemp, zsurf 
    72       REAL(wp) ::   zrhoa  = 1.22         ! Air density kg/m3 
    73       REAL(wp) ::   zcdrag = 1.5e-3       ! drag coefficient 
    74       REAL(wp) ::   ztx, zty, zmod, zcoef ! temporary variables 
    75       REAL(wp) ::   zyydd                 ! number of days in one year 
     68 
     69      INTEGER  ::   ios      ! Local integer output status for namelist read  
     70      !  
     71      CHARACTER(len=100) ::  cn_dir   !   Root directory for location of flux files  
     72      TYPE(FLD_N) ::   sn_wndi, sn_wndj                        ! informations about the fields to be read  
     73 
     74      NAMELIST/namsbc_usr/ ln_use_sbc, cn_dir , rn_vfac,  &  
     75         &                   sn_wndi, sn_wndj, rn_charn_const 
    7676      !!--------------------------------------------------------------------- 
    77       zyydd = REAL(nyear_len(1),wp) 
     77      !  
     78      IF( kt == nit000 ) THEN  
     79           
     80           
     81         REWIND( numnam_cfg )              ! Namelist namsbc_usr in configuration namelist  
     82         READ  ( numnam_cfg, namsbc_usr, IOSTAT = ios, ERR = 902 )  
     83902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_surge in configuration namelist' ) 
    7884 
    79       ! ---------------------------- ! 
    80       !  heat and freshwater fluxes  ! 
    81       ! ---------------------------- ! 
    82       !same temperature, E-P as in HAZELEGER 2000 
     85         IF(lwm) WRITE( numond, namsbc_usr ) 
     86         IF(lwp) WRITE(numout,*)' usr_sbc : AMM7_SURGE tide only case: NO surface forcing'  
     87         IF(lwp) WRITE(numout,*)' ~~~~~~~~~~~   utau = vtau = taum = wndm = qns = qsr = emp = sfx = 0' 
    8388 
    84       zyear0     =   ndate0 / 10000                             ! initial year 
    85       zmonth0    = ( ndate0 - zyear0 * 10000 ) / 100            ! initial month 
    86       zday0      =   ndate0 - zyear0 * 10000 - zmonth0 * 100    ! initial day betwen 1 and 30 
    87       zday_year0 = ( zmonth0 - 1 ) * 30.+zday0                  ! initial day betwen 1 and 360 
    88  
    89       ! current day (in hours) since january the 1st of the current year 
    90       ztime = REAL( kt ) * rdt / (rmmss * rhhmm)   &       !  total incrementation (in hours) 
    91          &      - (nyear  - 1) * rjjhh * zyydd             !  minus years since beginning of experiment (in hours) 
    92  
    93       ztimemax1 = ((5.*30.)+21.)* 24.                      ! 21th june     at 24h in hours 
    94       ztimemin1 = ztimemax1 + rjjhh * zyydd / 2            ! 21th december        in hours 
    95       ztimemax2 = ((6.*30.)+21.)* 24.                      ! 21th july     at 24h in hours 
    96       ztimemin2 = ztimemax2 - rjjhh * zyydd / 2            ! 21th january         in hours 
    97       !                                                    ! NB: rjjhh * zyydd / 4 = one seasonal cycle in hours 
    98  
    99       ! amplitudes 
    100       zemp_S    = 0.7       ! intensity of COS in the South 
    101       zemp_N    = 0.8       ! intensity of COS in the North 
    102       zemp_sais = 0.1 
    103       zTstar    = 28.3      ! intemsity from 28.3 a -5 deg 
    104  
    105       ! 1/2 period between 21th June and 21th December and between 21th July and 21th January 
    106       zcos_sais1 = COS( (ztime - ztimemax1) / (ztimemin1 - ztimemax1) * rpi )  
    107       zcos_sais2 = COS( (ztime - ztimemax2) / (ztimemax2 - ztimemin2) * rpi ) 
    108  
    109       ztrp= - 40.e0        ! retroaction term on heat fluxes (W/m2/K) 
    110       zconv = 3.16e-5      ! convertion factor: 1 m/yr => 3.16e-5 mm/s 
    111       DO jj = 1, jpj 
    112          DO ji = 1, jpi 
    113             ! domain from 15 deg to 50 deg between 27 and 28  degC at 15N, -3 
    114             ! and 13 degC at 50N 53.5 + or - 11 = 1/4 period : 
    115             ! 64.5 in summer, 42.5 in winter 
    116             t_star = zTstar * ( 1. + 1. / 50. * zcos_sais2 )                & 
    117                &                    * COS( rpi * (gphit(ji,jj) - 5.)               & 
    118                &                    / ( 53.5 * ( 1 + 11 / 53.5 * zcos_sais2 ) * 2.) ) 
    119             ! 23.5 deg : tropics 
    120             qsr (ji,jj) =  230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais1 ) / ( 0.9 * 180 ) ) 
    121             qns (ji,jj) = ztrp * ( tsb(ji,jj,1,jp_tem) - t_star ) - qsr(ji,jj) 
    122             IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj) ) THEN    ! zero at 37.8 deg, max at 24.6 deg 
    123                emp  (ji,jj) =   zemp_S * zconv   & 
    124                   &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) )  & 
    125                   &         * ( 1 - zemp_sais / zemp_S * zcos_sais1) 
    126             ELSE 
    127                emp (ji,jj) =  - zemp_N * zconv   & 
    128                   &         * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (46.8 - 37.2) )  & 
    129                   &         * ( 1 - zemp_sais / zemp_N * zcos_sais1 ) 
    130             ENDIF 
    131          END DO 
    132       END DO 
    133  
    134       zsumemp = GLOB_SUM( 'usrdef_sbc', emp  (:,:)   )  
    135       zsurf   = GLOB_SUM( 'usrdef_sbc', tmask(:,:,1) )  
    136       zsumemp = zsumemp / zsurf         ! Default GYRE configuration 
    137  
    138       ! freshwater (mass flux) and update of qns with heat content of emp 
    139       emp (:,:) = emp(:,:) - zsumemp * tmask(:,:,1)        ! freshwater flux (=0 in domain average) 
    140       sfx (:,:) = 0.0_wp                                   ! no salt flux 
    141       qns (:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp   ! evap and precip are at SST 
    142  
    143  
    144       ! ---------------------------- ! 
    145       !       momentum fluxes        ! 
    146       ! ---------------------------- ! 
    147       ! same wind as in Wico 
    148       !test date0 : ndate0 = 010203 
    149       zyear0  =   ndate0 / 10000 
    150       zmonth0 = ( ndate0 - zyear0 * 10000 ) / 100 
    151       zday0   =   ndate0 - zyear0 * 10000 - zmonth0 * 100 
    152       !Calculates nday_year, day since january 1st 
    153       zday_year0 = (zmonth0-1)*30.+zday0 
    154  
    155       !accumulates days of previous months of this year 
    156       ! day (in hours) since january the 1st 
    157       ztime = FLOAT( kt ) * rdt / (rmmss * rhhmm)  &  ! incrementation in hour 
    158          &     - (nyear - 1) * rjjhh * zyydd          !  - nber of hours the precedent years 
    159       ztimemax = ((5.*30.)+21.)* 24.               ! 21th june     in hours 
    160       ztimemin = ztimemax + rjjhh * zyydd / 2      ! 21th december in hours 
    161       !                                            ! NB: rjjhh * zyydd / 4 = 1 seasonal cycle in hours 
    162  
    163       ! mean intensity at 0.105 ; srqt(2) because projected with 45deg angle 
    164       ztau = 0.105 / SQRT( 2. ) 
    165       ! seasonal oscillation intensity 
    166       ztau_sais = 0.015 
    167       ztaun = ztau - ztau_sais * COS( (ztime - ztimemax) / (ztimemin - ztimemax) * rpi ) 
    168       DO jj = 1, jpj 
    169          DO ji = 1, jpi 
    170            ! domain from 15deg to 50deg and 1/2 period along 14deg 
    171            ! so 5/4 of half period with seasonal cycle 
    172            utau(ji,jj) = - ztaun * SIN( rpi * (gphiu(ji,jj) - 15.) / (29.-15.) ) 
    173            vtau(ji,jj) =   ztaun * SIN( rpi * (gphiv(ji,jj) - 15.) / (29.-15.) ) 
    174          END DO 
    175       END DO 
    176  
    177       ! module of wind stress and wind speed at T-point 
    178       zcoef = 1. / ( zrhoa * zcdrag )  
    179       DO jj = 2, jpjm1 
    180          DO ji = fs_2, fs_jpim1   ! vect. opt. 
    181             ztx = utau(ji-1,jj  ) + utau(ji,jj)  
    182             zty = vtau(ji  ,jj-1) + vtau(ji,jj)  
    183             zmod = 0.5 * SQRT( ztx * ztx + zty * zty ) 
    184             taum(ji,jj) = zmod 
    185             wndm(ji,jj) = SQRT( zmod * zcoef ) 
    186          END DO 
    187       END DO 
    188       CALL lbc_lnk_multi( 'usrdef_sbc', taum(:,:), 'T', 1. , wndm(:,:), 'T', 1. ) 
    189  
    190       ! ---------------------------------- ! 
    191       !  control print at first time-step  ! 
    192       ! ---------------------------------- ! 
    193       IF( kt == nit000 .AND. lwp ) THEN  
    194          WRITE(numout,*) 
    195          WRITE(numout,*)'usrdef_sbc_oce : analytical surface fluxes for GYRE configuration'                
    196          WRITE(numout,*)'~~~~~~~~~~~ '  
    197          WRITE(numout,*)'           nyear      = ', nyear 
    198          WRITE(numout,*)'           nmonth     = ', nmonth 
    199          WRITE(numout,*)'           nday       = ', nday 
    200          WRITE(numout,*)'           nday_year  = ', nday_year 
    201          WRITE(numout,*)'           ztime      = ', ztime 
    202          WRITE(numout,*)'           ztimemax   = ', ztimemax 
    203          WRITE(numout,*)'           ztimemin   = ', ztimemin 
    204          WRITE(numout,*)'           ztimemax1  = ', ztimemax1 
    205          WRITE(numout,*)'           ztimemin1  = ', ztimemin1 
    206          WRITE(numout,*)'           ztimemax2  = ', ztimemax2 
    207          WRITE(numout,*)'           ztimemin2  = ', ztimemin2 
    208          WRITE(numout,*)'           zyear0     = ', zyear0 
    209          WRITE(numout,*)'           zmonth0    = ', zmonth0 
    210          WRITE(numout,*)'           zday0      = ', zday0 
    211          WRITE(numout,*)'           zday_year0 = ', zday_year0 
    212          WRITE(numout,*)'           zyydd      = ', zyydd 
    213          WRITE(numout,*)'           zemp_S     = ', zemp_S 
    214          WRITE(numout,*)'           zemp_N     = ', zemp_N 
    215          WRITE(numout,*)'           zemp_sais  = ', zemp_sais 
    216          WRITE(numout,*)'           zTstar     = ', zTstar 
    217          WRITE(numout,*)'           zsumemp    = ', zsumemp 
    218          WRITE(numout,*)'           zsurf      = ', zsurf 
    219          WRITE(numout,*)'           ztrp       = ', ztrp 
    220          WRITE(numout,*)'           zconv      = ', zconv 
    221          WRITE(numout,*)'           ndastp     = ', ndastp 
    222          WRITE(numout,*)'           adatrj     = ', adatrj 
     89         utau(:,:) = 0._wp  
     90         vtau(:,:) = 0._wp  
     91         taum(:,:) = 0._wp  
     92         wndm(:,:) = 0._wp  
     93         !  
     94         emp (:,:) = 0._wp  
     95         sfx (:,:) = 0._wp  
     96         qns (:,:) = 0._wp  
     97         qsr (:,:) = 0._wp  
     98         ! 
    22399      ENDIF 
    224100      ! 
     
    231107 
    232108 
    233    SUBROUTINE usrdef_sbc_ice_flx( kt, phs, phi ) 
     109   SUBROUTINE usrdef_sbc_ice_flx( kt ) 
    234110      INTEGER, INTENT(in) ::   kt   ! ocean time step 
    235       REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness 
    236       REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness 
    237111   END SUBROUTINE usrdef_sbc_ice_flx 
    238112 
  • NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/par_oce.F90

    r14075 r15422  
    6464   INTEGER, PUBLIC, PARAMETER ::   jp_tem = 1    !: indice for temperature 
    6565   INTEGER, PUBLIC, PARAMETER ::   jp_sal = 2    !: indice for salinity 
     66   INTEGER, PUBLIC, PARAMETER ::   jp_dep = 3    !: indice for depth   
     67   INTEGER, PUBLIC, PARAMETER ::   jp_msk = 4    !: indice for mask 
    6668 
    6769   !!---------------------------------------------------------------------- 
  • NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/step.F90

    r14075 r15422  
    206206                         CALL dia_ar5 ( kstp )        ! ar5 diag 
    207207      IF( ln_diaptr  )   CALL dia_ptr                 ! Poleward adv/ldf TRansports diagnostics 
     208      IF( lk_diaharm_fast )                           &   
     209            &            CALL dia_harm_fast( kstp )   ! Tidal harmonic analysis - restart and faster version  
    208210      IF( ln_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis 
    209211                         CALL dia_wri ( kstp )        ! ocean model: outputs 
  • NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/step_oce.F90

    r14075 r15422  
    7979   USE diahth          ! thermocline depth                (dia_hth routine) 
    8080   USE diahsb          ! heat, salt and volume budgets    (dia_hsb routine) 
     81   USE diaharm_fast    ! harmonic analysis of tides (harm_ana routine) 
    8182   USE diaharm 
    8283   USE diacfl 
  • NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/stpctl.F90

    r14075 r15422  
    126126         zmax(1) = MAXVAL(  ABS( sshn(:,:) )  )                               ! ssh max 
    127127      ENDIF 
    128       zmax(2) = MAXVAL(  ABS( un(:,:,:) )  )                                  ! velocity max (zonal only) 
     128      zmax(2) = MAXVAL(  un(:,:,:)*un(:,:,:) + vn(:,:,:)*vn(:,:,:)  )         ! velocity max 
    129129      llmsk(:,:,:) = tmask(:,:,:) == 1._wp 
    130130      IF( COUNT( llmsk(:,:,:) ) > 0 ) THEN   ! avoid huge values sent back for land processors...       
     
    167167      !                                   !==  error handling  ==! 
    168168      IF(   zmax(1) >   20._wp .OR.   &                    ! too large sea surface height ( > 20 m ) 
    169          &  zmax(2) >   10._wp .OR.   &                    ! too large velocity ( > 10 m/s) 
     169         &  zmax(2) >  100._wp .OR.   &                    ! too large velocity ( > 10 m/s) 
    170170         &  zmax(3) >=   0._wp .OR.   &                    ! negative or zero sea surface salinity 
    171171         &  zmax(4) >= 100._wp .OR.   &                    ! too large sea surface salinity ( > 100 ) 
     
    177177            IF( lwm .AND. kt /= nitend )   istatus = NF90_CLOSE(idrun) 
    178178            CALL mpp_maxloc( 'stpctl', ABS(sshn)        , ssmask(:,:)  , zzz, ih(1:2)  )   ;   ih(3) = 0 
    179             CALL mpp_maxloc( 'stpctl', ABS(un)          , umask (:,:,:), zzz, iu  ) 
     179            CALL mpp_maxloc( 'stpctl', un*un + vn*vn    , umask (:,:,:), zzz, iu  ) 
    180180            CALL mpp_minloc( 'stpctl', tsn(:,:,:,jp_sal), tmask (:,:,:), zzz, is1 ) 
    181181            CALL mpp_maxloc( 'stpctl', tsn(:,:,:,jp_sal), tmask (:,:,:), zzz, is2 ) 
     
    192192         ELSE 
    193193            ih(1:2)= MAXLOC( ABS( sshn(:,:)   )                              ) + (/ nimpp - 1, njmpp - 1    /)   ;   ih(3) = 0 
    194             iu(:)  = MAXLOC( ABS( un  (:,:,:) )                              ) + (/ nimpp - 1, njmpp - 1, 0 /) 
     194            iu(:)  = MAXLOC( un(:,:,:)*un(:,:,:) + vn(:,:,:)*vn(:,:,:)       ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    195195            is1(:) = MINLOC( tsn(:,:,:,jp_sal), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
    196196            is2(:) = MAXLOC( tsn(:,:,:,jp_sal), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 
     
    200200         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests' 
    201201         CALL wrt_line(ctmp2, kt, ' |ssh| max ',   zmax(1), ih , iareasum(1), iareamin(1), iareamax(1) )  
    202          CALL wrt_line(ctmp3, kt, ' |U|   max ',   zmax(2), iu , iareasum(2), iareamin(2), iareamax(2) )  
     202         CALL wrt_line(ctmp3, kt, ' Vel2  max ',   zmax(2), iu , iareasum(2), iareamin(2), iareamax(2) )  
    203203         CALL wrt_line(ctmp4, kt, ' Sal   min ', - zmax(3), is1, iareasum(3), iareamin(3), iareamax(3) )  
    204204         CALL wrt_line(ctmp5, kt, ' Sal   max ',   zmax(4), is2, iareasum(4), iareamin(4), iareamax(4) )  
     
    226226      ENDIF 
    227227      ! 
    228 9500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16) 
     2289500  FORMAT(' it :', i8, '    |ssh|_max: ', D23.16, ' Vel2_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16) 
    229229      ! 
    230230   END SUBROUTINE stp_ctl 
Note: See TracChangeset for help on using the changeset viewer.