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 901 for trunk/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90 – NEMO

Ignore:
Timestamp:
2008-04-22T21:16:18+02:00 (16 years ago)
Author:
rblod
Message:

Change ldf computation to take into account the grid anisotropy, see ticket #115

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OPA_SRC/LDF/ldfdyn_c3d.h90

    r719 r901  
    3838         !               ! ( 0 < zr < 1 ) 
    3939         zh = 500.,   &  ! depth of at which start the reduction ( > dept(1) ) 
    40          zdx_max  ,   &  ! maximum grid spacing over the global domain 
     40         zd_max   ,   &  ! maximum grid spacing over the global domain 
    4141         za00, zc, zd    ! temporary scalars 
     42      REAL(wp) ::        & 
     43         zetmax, zefmax, & 
     44         zeumax, zevmax    
    4245      REAL(wp), DIMENSION(jpk) ::   zcoef   ! temporary workspace 
    4346      !!---------------------------------------------------------------------- 
     
    5457         ! (USER: modify ahm1 and ahm2 following your desiderata) 
    5558 
    56          zdx_max = MAXVAL( e1t(:,:) ) 
    57          IF( lk_mpp )   CALL mpp_max( zdx_max )   ! max over the global domain 
     59         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
     60         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain 
    5861 
    5962         IF(lwp) WRITE(numout,*) '              laplacian operator: ahm proportional to e1' 
    60          IF(lwp) WRITE(numout,*) '              Caution, here we assume your mesh is isotropic ...' 
    61          IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zdx_max, ' maximum value for ahm = ', ahm0 
    62  
    63  
    64          za00 = ahm0 / zdx_max 
     63         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0 
     64 
     65         za00 = ahm0 / zd_max 
    6566 
    6667         IF( ln_dynldf_iso ) THEN 
     
    7475         CALL ldf_zpf( .TRUE. , 1000., 500., 0.25, fsdept(:,:,:), ahm2 )   ! vertical profile 
    7576         DO jk = 1,jpk 
    76             ahm1(:,:,jk) = za00 * e1t(:,:) * ahm1(:,:,jk) 
    77             ahm2(:,:,jk) = za00 * e1f(:,:) * ahm2(:,:,jk) 
     77            DO jj = 1, jpj 
     78               DO ji = 1, jpi 
     79                  zetmax = MAX( e1t(ji,jj), e2t(ji,jj) ) 
     80                  zefmax = MAX( e1f(ji,jj), e2f(ji,jj) ) 
     81                  ahm1(ji,jj,jk) = za00 * zetmax * ahm1(ji,jj,jk) 
     82                  ahm2(ji,jj,jk) = za00 * zefmax * ahm2(ji,jj,jk) 
     83               END DO 
     84            END DO 
    7885         END DO 
    7986 
     
    109116      IF( ln_dynldf_bilap ) THEN 
    110117 
    111          zdx_max = MAXVAL( e1u(:,:) ) 
    112          IF( lk_mpp )   CALL mpp_max( zdx_max )   ! max over the global domain 
     118         zd_max = MAX( MAXVAL( e1u(:,:) ), MAXVAL( e2u(:,:) ) ) 
     119         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain 
    113120 
    114121         IF(lwp) WRITE(numout,*) '              bi-laplacian operator: ahm proportional to e1**3 ' 
    115          IF(lwp) WRITE(numout,*) '              Caution, here we assume your mesh is isotropic ...' 
    116          IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zdx_max, ' maximum value for ahm = ', ahm0 
    117  
    118          za00 = ahm0 / ( zdx_max * zdx_max * zdx_max ) 
    119          ahm3(:,:,1) = za00 * e1u(:,:) * e1u(:,:) * e1u(:,:) 
    120          ahm4(:,:,1) = za00 * e1v(:,:) * e1v(:,:) * e1v(:,:) 
     122         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0 
     123 
     124         za00 = ahm0 / ( zd_max * zd_max * zd_max ) 
     125         DO jj = 1, jpj 
     126            DO ji = 1, jpi 
     127               zeumax = MAX( e1u(ji,jj), e2u(ji,jj) ) 
     128               zevmax = MAX( e1v(ji,jj), e2v(ji,jj) ) 
     129               ahm3(ji,jj,1) = za00 * zeumax * zeumax * zeumax 
     130               ahm4(ji,jj,1) = za00 * zevmax * zevmax * zevmax 
     131            END DO 
     132         END DO 
    121133 
    122134         zh = MAX( zh, fsdept(1,1,1) )   ! at least the first reach ahm0 
Note: See TracChangeset for help on using the changeset viewer.