Changeset 901


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

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

Location:
trunk/NEMO/OPA_SRC/LDF
Files:
4 edited

Legend:

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

    r719 r901  
    3636 
    3737      !! * Local variables 
    38       REAL(wp) ::   za00, zdx_max 
     38      INTEGER :: ji, jj 
     39      REAL(wp) ::   za00, zd_max, zetmax, zeumax, zefmax, zevmax 
    3940      !!---------------------------------------------------------------------- 
    4041 
     
    5051         ! (USER: modify ahm1 and ahm2 following your desiderata) 
    5152 
    52          zdx_max = MAXVAL( e1t(:,:) ) 
    53          IF( lk_mpp )   CALL mpp_max( zdx_max )   ! max over the global domain 
     53         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
     54         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain 
    5455 
    5556         IF(lwp) WRITE(numout,*) '              laplacian operator: ahm proportional to e1' 
    56          IF(lwp) WRITE(numout,*) '              Caution, here we assume your mesh is isotropic ...' 
    57          IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zdx_max, ' maximum value for ahm = ', ahm0 
    58  
    59          za00 = ahm0 / zdx_max 
    60          ahm1(:,:) = za00 * e1t(:,:) 
    61          ahm2(:,:) = za00 * e1f(:,:) 
     57         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0 
     58 
     59         za00 = ahm0 / zd_max 
     60         DO jj = 1, jpj 
     61            DO ji = 1, jpi 
     62               zetmax = MAX( e1t(ji,jj), e2t(ji,jj) ) 
     63               zefmax = MAX( e1f(ji,jj), e2f(ji,jj) ) 
     64               ahm1(ji,jj) = za00 * zetmax 
     65               ahm2(ji,jj) = za00 * zefmax 
     66            END DO 
     67         END DO 
    6268 
    6369         IF( ln_dynldf_iso ) THEN 
     
    9096         !       in the to horizontal direction 
    9197 
    92          zdx_max = MAXVAL( e1u(:,:) ) 
    93          IF( lk_mpp )   CALL mpp_max( zdx_max )   ! max over the global domain 
     98         zd_max = MAX( MAXVAL( e1u(:,:) ), MAXVAL( e2u(:,:) ) ) 
     99         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain 
    94100 
    95101         IF(lwp) WRITE(numout,*) '              bi-laplacian operator: ahm proportional to e1**3 ' 
    96          IF(lwp) WRITE(numout,*) '              Caution, here we assume your mesh is isotropic ...' 
    97          IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zdx_max, ' maximum value for ahm = ', ahm0 
    98  
    99          za00 = ahm0 / ( zdx_max * zdx_max * zdx_max ) 
    100          ahm3(:,:) = za00 * e1u(:,:) * e1u(:,:) * e1u(:,:) 
    101          ahm4(:,:) = za00 * e1v(:,:) * e1v(:,:) * e1v(:,:) 
     102         IF(lwp) WRITE(numout,*) '              maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0 
     103 
     104         za00 = ahm0 / ( zd_max * zd_max * zd_max ) 
     105         DO jj = 1, jpj 
     106            DO ji = 1, jpi 
     107               zeumax = MAX( e1u(ji,jj), e2u(ji,jj) ) 
     108               zevmax = MAX( e1v(ji,jj), e2v(ji,jj) ) 
     109               ahm3(ji,jj) = za00 * zeumax * zeumax * zeumax 
     110               ahm4(ji,jj) = za00 * zevmax * zevmax * zevmax 
     111            END DO 
     112         END DO 
    102113 
    103114         ! Control print 
  • 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 
  • trunk/NEMO/OPA_SRC/LDF/ldftra_c2d.h90

    r719 r901  
    3131 
    3232      !! * Local variables 
    33 # if defined key_traldf_eiv && defined key_orca_r4 
    3433      INTEGER ::   ji, jj                  ! dummy loop indices 
    35 # endif 
    3634# if defined key_orca_r4 
    3735      INTEGER :: i1, i2, j1, j2 
    3836# endif 
    39       REAL(wp) ::   za00, zdx_max 
     37      REAL(wp) ::   za00, zd_max, zeumax, zevmax, zetmax 
    4038       
    4139      !!---------------------------------------------------------------------- 
     
    5351      ENDIF 
    5452 
    55       zdx_max = MAXVAL( e1t(:,:) ) 
    56       IF( lk_mpp )   CALL mpp_max( zdx_max )   ! max over the global domain 
     53      zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
     54      IF( lk_mpp ) CALL mpp_max( zd_max )   ! max over the global domain 
    5755 
    5856      ! harmonic operator : (U-, V-, W-points) 
     
    6058      IF( ln_traldf_lap ) THEN 
    6159   
    62          za00 = aht0 / zdx_max 
     60         za00 = aht0 / zd_max 
    6361   
    64 !RB         ahtu(:,:) = za00 * e1u(:,:)               ! set ahtu = ahtv at u- and v-points, 
    65 !RB         ahtv(:,:) = za00 * e1f(:,:)               ! and ahtw at w-point (idem T-point) 
    66 !RB         ahtw(:,:) = za00 * e1t(:,:)               !  
    67          ahtu(:,:) = aht0                ! set ahtu = ahtv at u- and v-points, 
    68          ahtv(:,:) = aht0                ! and ahtw at w-point (idem T-point) 
    69          ahtw(:,:) = aht0                ! 
    70  
    71  
     62         DO jj = 1, jpj  
     63            DO ji = 1, jpi  
     64               zeumax = MAX( e1u(ji,jj), e2u(ji,jj) )  
     65               zevmax = MAX( e1v(ji,jj), e2v(ji,jj) )  
     66               zetmax = MAX( e1t(ji,jj), e2t(ji,jj) ) 
     67               ahtu(ji,jj) = za00 * zeumax ! set ahtu = ahtv at u- and v-points,  
     68               ahtv(ji,jj) = za00 * zevmax ! and ahtw at w-point (idem T-point)  
     69               ahtw(ji,jj) = za00 * zetmax !  
     70            END DO 
     71         END DO 
    7272 
    7373         CALL lbc_lnk( ahtu, 'U', 1. )   ! Lateral boundary conditions 
    7474         CALL lbc_lnk( ahtv, 'V', 1. )   ! (no change of sign) 
    7575         CALL lbc_lnk( ahtw, 'W', 1. ) 
     76 
     77         ! Special case for ORCA R2 and R4 configurations (overwrite the value of ahtu ahtv ahtw) 
     78         ! ============================================== 
     79         IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) )   THEN 
     80            ahtu(:,:) = aht0              ! set ahtu = ahtv at u- and v-points, 
     81            ahtv(:,:) = aht0              ! and ahtw at w-point 
     82            ahtw(:,:) = aht0              ! (here : no space variation) 
     83            IF(lwp) WRITE(numout,*) '               ORCA R2 or R4 case' 
     84            IF(lwp) WRITE(numout,*) '               Constant values used for eddy diffusivity coefficients' 
     85            IF(lwp) WRITE(numout,*) '               Variation lat/lon only for eddy induced velocity coefficients' 
     86         ENDIF 
    7687 
    7788         ! Control print 
     
    99110         !       in the to horizontal direction 
    100111 
    101          za00 = aht0 / ( zdx_max * zdx_max * zdx_max ) 
    102          ahtt(:,:) = za00 * e1t(:,:) * e1t(:,:) *e1t(:,:)      ! set ahtt at T-point  
     112         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 
     113         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain 
     114 
     115         za00 = aht0 / ( zd_max * zd_max * zd_max ) 
     116         DO jj = 1, jpj 
     117            DO ji = 1, jpi 
     118               zetmax = MAX( e1t(ji,jj), e2t(ji,jj) ) 
     119               ahtt(ji,jj) = za00 * zetmax * zetmax * zetmax      ! set ahtt at T-point 
     120            END DO 
     121         END DO 
    103122 
    104123         CALL lbc_lnk( ahtt, 'T', 1. )   ! Lateral boundary conditions on ( ahtt ) 
  • trunk/NEMO/OPA_SRC/LDF/ldftra_c3d.h90

    r719 r901  
    4545         IF(lwp) WRITE(numout,*) ' ldf_tra_c3d : 3D eddy diffusivity and eddy induced velocity coefficients' 
    4646         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~   --  ' 
     47         IF(lwp) WRITE(numout,*) '               Coefficients set to constant' 
    4748         IF(lwp) WRITE(numout,*) 
    4849      ELSE 
     
    5051         IF(lwp) WRITE(numout,*) ' ldf_tra_c3d : 3D eddy diffusivity coefficient' 
    5152         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~   -- ' 
     53         IF(lwp) WRITE(numout,*) '               Coefficients set to constant' 
    5254         IF(lwp) WRITE(numout,*) 
    5355      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.