!!---------------------------------------------------------------------- !! *** ldfdyn_c3d.h90 *** !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! OPA 9.0 , LOCEAN-IPSL (2005) !! $Id$ !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! 'key_dynldf_c3d' 3D lateral eddy viscosity coefficients !!---------------------------------------------------------------------- SUBROUTINE ldf_dyn_c3d( ld_print ) !!---------------------------------------------------------------------- !! *** ROUTINE ldf_dyn_c3d *** !! !! ** Purpose : initializations of the horizontal ocean physics !! !! ** Method : 3D eddy viscosity coef. ( longitude, latitude, depth ) !! laplacian operator : ahm1, ahm2 defined at T- and F-points !! ahm2, ahm4 never used !! bilaplacian operator : ahm1, ahm2 never used !! : ahm3, ahm4 defined at U- and V-points !! ??? explanation of the default is missing !!---------------------------------------------------------------------- !! * Modules used USE ldftra_oce, ONLY : aht0 !! * Arguments LOGICAL, INTENT (in) :: ld_print ! If true, output arrays on numout !! * local variables INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: & zr = 0.2 , & ! maximum of the reduction factor at the bottom ocean ! ! ( 0 < zr < 1 ) zh = 500., & ! depth of at which start the reduction ( > dept(1) ) zd_max , & ! maximum grid spacing over the global domain za00, zc, zd ! temporary scalars REAL(wp) :: & zetmax, zefmax, & zeumax, zevmax REAL(wp), DIMENSION(jpk) :: zcoef ! temporary workspace !!---------------------------------------------------------------------- IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'ldf_dyn_c3d : 3D lateral eddy viscosity coefficient' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' ! Set ahm1 and ahm2 ( T- and F- points) (used for laplacian operators ! ================= whatever its orientation is) IF( ln_dynldf_lap ) THEN ! define ahm1 and ahm2 at the right grid point position ! (USER: modify ahm1 and ahm2 following your desiderata) zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain IF(lwp) WRITE(numout,*) ' laplacian operator: ahm proportional to e1' IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0 za00 = ahm0 / zd_max IF( ln_dynldf_iso ) THEN IF(lwp) WRITE(numout,*) ' Caution, as implemented now, the isopycnal part of momentum' IF(lwp) WRITE(numout,*) ' mixing use aht0 as eddy viscosity coefficient. Thus, it is' IF(lwp) WRITE(numout,*) ' uniform and you must be sure that your ahm is greater than' IF(lwp) WRITE(numout,*) ' aht0 everywhere in the model domain.' ENDIF CALL ldf_zpf( .TRUE. , 1000., 500., 0.25, fsdept(:,:,:), ahm1 ) ! vertical profile CALL ldf_zpf( .TRUE. , 1000., 500., 0.25, fsdept(:,:,:), ahm2 ) ! vertical profile DO jk = 1,jpk DO jj = 1, jpj DO ji = 1, jpi zetmax = MAX( e1t(ji,jj), e2t(ji,jj) ) zefmax = MAX( e1f(ji,jj), e2f(ji,jj) ) ahm1(ji,jj,jk) = za00 * zetmax * ahm1(ji,jj,jk) ahm2(ji,jj,jk) = za00 * zefmax * ahm2(ji,jj,jk) END DO END DO END DO ! Special case for ORCA R2 and R4 configurations (overwrite the value of ahm1 ahm2) ! ============================================== IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' ORCA R2 or R4: overwrite the previous definition of ahm' IF(lwp) WRITE(numout,*) ' =============' CALL ldf_dyn_c3d_orca( ld_print ) ENDIF ENDIF ! Control print IF(lwp .AND. ld_print ) THEN WRITE(numout,*) WRITE(numout,*) ' 3D ahm1 array (k=1)' CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout ) WRITE(numout,*) WRITE(numout,*) ' 3D ahm2 array (k=1)' CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout ) ENDIF ! ahm3 and ahm4 at U- and V-points (used for bilaplacian operator ! ================================ whatever its orientation is) ! (USER: modify ahm3 and ahm4 following your desiderata) ! Here: ahm is proportional to the cube of the maximum of the gridspacing ! in the to horizontal direction IF( ln_dynldf_bilap ) THEN zd_max = MAX( MAXVAL( e1u(:,:) ), MAXVAL( e2u(:,:) ) ) IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain IF(lwp) WRITE(numout,*) ' bi-laplacian operator: ahm proportional to e1**3 ' IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahm = ', ahm0 za00 = ahm0_blp / ( zd_max * zd_max * zd_max ) DO jj = 1, jpj DO ji = 1, jpi zeumax = MAX( e1u(ji,jj), e2u(ji,jj) ) zevmax = MAX( e1v(ji,jj), e2v(ji,jj) ) ahm3(ji,jj,1) = za00 * zeumax * zeumax * zeumax ahm4(ji,jj,1) = za00 * zevmax * zevmax * zevmax END DO END DO zh = MAX( zh, fsdept(1,1,1) ) ! at least the first reach ahm0 IF( ln_zco ) THEN ! z-coordinate, same profile everywhere IF(lwp) WRITE(numout,'(36x," ahm ", 7x)') DO jk = 1, jpk IF( fsdept(1,1,jk) <= zh ) THEN zcoef(jk) = 1.e0 ELSE zcoef(jk) = 1.e0 + ( zr - 1.e0 ) & & * ( 1. - EXP( ( fsdept(1,1,jk ) - zh ) / zh ) ) & & / ( 1. - EXP( ( fsdept(1,1,jpkm1) - zh ) / zh ) ) ENDIF ahm3(:,:,jk) = ahm3(:,:,1) * zcoef(jk) ahm4(:,:,jk) = ahm4(:,:,1) * zcoef(jk) IF(lwp) WRITE(numout,'(34x,E7.2,8x,i3)') zcoef(jk) * ahm0, jk END DO ELSE ! partial steps or s-ccordinate # if defined key_zco zc = gdept_0(jpkm1) # else zc = MAXVAL( fsdept(:,:,jpkm1) ) # endif IF( lk_mpp ) CALL mpp_max( zc ) ! max over the global domain zc = 1. / ( 1. - EXP( ( zc - zh ) / zh ) ) DO jk = 2, jpkm1 DO jj = 1, jpj DO ji = 1, jpi IF( fsdept(ji,jj,jk) <= zh ) THEN ahm3(ji,jj,jk) = ahm3(ji,jj,1) ahm4(ji,jj,jk) = ahm4(ji,jj,1) ELSE zd = 1.e0 + ( zr - 1.e0 ) * ( 1. - EXP( ( fsdept(ji,jj,jk) - zh ) / zh ) ) * zc ahm3(ji,jj,jk) = ahm3(ji,jj,1) * zd ahm4(ji,jj,jk) = ahm4(ji,jj,1) * zd ENDIF END DO END DO END DO ahm3(:,:,jpk) = ahm3(:,:,jpkm1) ahm4(:,:,jpk) = ahm4(:,:,jpkm1) IF(lwp) WRITE(numout,'(36x," ahm ", 7x)') DO jk = 1, jpk IF(lwp) WRITE(numout,'(30x,E10.2,8x,i3)') ahm3(1,1,jk), jk END DO ENDIF ! Control print IF( lwp .AND. ld_print ) THEN WRITE(numout,*) WRITE(numout,*) 'inildf: ahm3 array at level 1' CALL prihre(ahm3(:,:,1 ),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) WRITE(numout,*) WRITE(numout,*) 'inildf: ahm4 array at level 1' CALL prihre(ahm4(:,:,jpk),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) ENDIF ENDIF END SUBROUTINE ldf_dyn_c3d SUBROUTINE ldf_dyn_c3d_orca( ld_print ) !!---------------------------------------------------------------------- !! *** ROUTINE ldf_dyn_c3d *** !! !! ** Purpose : ORCA R2 an R4 only !! !! ** Method : blah blah blah .... !!---------------------------------------------------------------------- !! * Modules used USE ldftra_oce, ONLY : aht0 !! * Arguments LOGICAL, INTENT (in) :: ld_print ! If true, output arrays on numout !! * local variables INTEGER :: ji, jj, jk, jn ! dummy loop indices INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers INTEGER :: inum ! temporary logical unit INTEGER :: iim, ijm INTEGER :: ifreq, il1, il2, ij, ii INTEGER, DIMENSION(jpidta, jpjdta) :: idata INTEGER, DIMENSION(jpi , jpj ) :: icof REAL(wp) :: & zahmeq, zcoff, zcoft, zmsk, & ! ??? zemax, zemin, zeref, zahmm REAL(wp), DIMENSION(jpi,jpj) :: zahm0 REAL(wp), DIMENSION(jpk) :: zcoef CHARACTER (len=15) :: clexp !!---------------------------------------------------------------------- IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'ldfdyn_c3d_orca : 3D eddy viscosity coefficient' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~' IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' orca R2 or R4 ocean model' IF(lwp) WRITE(numout,*) ' reduced in the surface Eq. strip ' IF(lwp) WRITE(numout,*) ! Read 2d integer array to specify western boundary increase in the ! ===================== equatorial strip (20N-20S) defined at t-points CALL ctl_opn( inum, 'ahmcoef', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) READ(inum,9101) clexp, iim, ijm READ(inum,'(/)') ifreq = 40 il1 = 1 DO jn = 1, jpidta/ifreq+1 READ(inum,'(/)') il2 = MIN( jpidta, il1+ifreq-1 ) READ(inum,9201) ( ii, ji = il1, il2, 5 ) READ(inum,'(/)') DO jj = jpjdta, 1, -1 READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 ) END DO il1 = il1 + ifreq END DO DO jj = 1, nlcj DO ji = 1, nlci icof(ji,jj) = idata( mig(ji), mjg(jj) ) END DO END DO DO jj = nlcj+1, jpj DO ji = 1, nlci icof(ji,jj) = icof(ji,nlcj) END DO END DO DO jj = 1, jpj DO ji = nlci+1, jpi icof(ji,jj) = icof(nlci,jj) END DO END DO 9101 FORMAT(1x,a15,2i8) 9201 FORMAT(3x,13(i3,12x)) 9202 FORMAT(i3,41i3) ! Set ahm1 and ahm2 ! ================= ! define ahm1 and ahm2 at the right grid point position ! (USER: modify ahm1 and ahm2 following your desiderata) ! biharmonic : ahm1 (ahm2) defined at u- (v-) point ! harmonic : ahm1 (ahm2) defined at t- (f-) point ! first level : as for 2D coefficients ! Decrease ahm to zahmeq m2/s in the tropics ! (from 90 to 20 degre: ahm = constant ! from 20 to 2.5 degre: ahm = decrease in (1-cos)/2 ! from 2.5 to 0 degre: ahm = constant ! symmetric in the south hemisphere) IF( jp_cfg == 4 ) THEN zahmeq = 5.0 * aht0 zahmm = min( 160000.0, ahm0) zemax = MAXVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 ) zemin = MINVAL ( e1t(:,:) * e2t(:,:), tmask(:,:,1) .GE. 0.5 ) zeref = MAXVAL ( e1t(:,:) * e2t(:,:), & & tmask(:,:,1) .GE. 0.5 .AND. ABS(gphit(:,:)) .GT. 50. ) DO jj = 1, jpj DO ji = 1, jpi zmsk = e1t(ji,jj) * e2t(ji,jj) IF( abs(gphit(ji,jj)) .LE. 15 ) THEN zahm0(ji,jj) = ahm0 ELSE IF( zmsk .GE. zeref ) THEN zahm0(ji,jj) = ahm0 ELSE zahm0(ji,jj) = zahmm + (ahm0-zahmm)*(1.0 - & & cos((rpi*0.5*(zmsk-zemin)/(zeref-zemin)))) ENDIF ENDIF END DO END DO ENDIF IF( jp_cfg == 2 ) THEN zahmeq = aht0 zahmm = ahm0 zahm0(:,:) = ahm0 ENDIF DO jj = 1, jpj DO ji = 1, jpi IF( ABS(gphif(ji,jj)) >= 20.) THEN ahm2(ji,jj,1) = zahm0(ji,jj) ELSEIF( ABS(gphif(ji,jj)) <= 2.5) THEN ahm2(ji,jj,1) = zahmeq ELSE ahm2(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2. & & *(1.-COS( rad*(ABS(gphif(ji,jj))-2.5)*180./17.5 ) ) ENDIF IF( ABS(gphit(ji,jj)) >= 20.) THEN ahm1(ji,jj,1) = zahm0(ji,jj) ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ahm1(ji,jj,1) = zahmeq ELSE ahm1(ji,jj,1) = zahmeq + (zahm0(ji,jj)-zahmeq)/2. & & *(1.-COS( rad*(ABS(gphit(ji,jj))-2.5)*180./17.5 ) ) ENDIF END DO END DO ! increase along western boundaries of equatorial strip ! t-point DO jj = 1, jpjm1 DO ji = 1, jpim1 zcoft = float( icof(ji,jj) ) / 100. ahm1(ji,jj,1) = zcoft * zahm0(ji,jj) + (1.-zcoft) * ahm1(ji,jj,1) END DO END DO ! f-point icof(:,:) = icof(:,:) * tmask(:,:,1) DO jj = 1, jpjm1 DO ji = 1, jpim1 ! NO vector opt. zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1) IF( zmsk == 0. ) THEN zcoff = 1. ELSE zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) ) & / (zmsk * 100.) ENDIF ahm2(ji,jj,1) = zcoff * zahm0(ji,jj) + (1.-zcoff) * ahm2(ji,jj,1) END DO END DO ! other level: re-increase the coef in the deep ocean #if defined key_orca_lev10 DO jk = 1, 210 zcoef(jk) = 1. END DO DO jk= 211, 230 zcoef(jk) = 1. + 0.1 * FLOAT(jk-210) END DO DO jk= 231, 260 zcoef(jk) = 3. + 0.2 * FLOAT(jk-230) END DO DO jk= 261, 270 zcoef(jk) = 9. + 0.1 * FLOAT(jk-260) END DO DO jk= 271, jpk zcoef(jk) = 10. END DO DO jk= 1, jpk IF(lwp) WRITE(numout,*) 'k= ',jk, 'cof ', zcoef(jk) END DO #else DO jk = 1, 21 zcoef(jk) = 1. END DO zcoef(22) = 2. zcoef(23) = 3. zcoef(24) = 5. zcoef(25) = 7. zcoef(26) = 9. DO jk = 27, jpk zcoef(jk) = 10. END DO #endif DO jk = 2, jpk ahm1(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm1(:,:,1) ) ahm2(:,:,jk) = MIN( zahm0(:,:), zcoef(jk) * ahm2(:,:,1) ) END DO IF( jp_cfg == 4 ) THEN ! Limit AHM in Gibraltar strait ij0 = 50 ; ij1 = 53 ii0 = 69 ; ii1 = 71 DO jk = 1, jpk ahm1(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk) = min( zahmm, ahm1 (mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk) ) ahm2(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk) = min( zahmm, ahm2 (mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , jk) ) END DO ENDIF ! Lateral boundary conditions on ( ahm1, ahm2 ) ! ============== CALL lbc_lnk( ahm1, 'T', 1. ) ! T-point, unchanged sign CALL lbc_lnk( ahm2, 'F', 1. ) ! F-point, unchanged sign ! Control print IF(lwp) THEN WRITE(numout,*) WRITE(numout,*) ' 3D ahm1 array (k=1)' CALL prihre( ahm1(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout ) WRITE(numout,*) WRITE(numout,*) ' 3D ahm2 array (k=1)' CALL prihre( ahm2(:,:,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout ) WRITE(numout,*) WRITE(numout,*) ' 3D ahm2 array (k=jpk)' CALL prihre( ahm2(:,:,jpk), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1.e-3, numout ) ENDIF ! Set ahm3 and ahm4 ! ================= ! define ahm3 and ahm4 at the right grid point position ! initialization to a constant value ! (USER: modify ahm3 and ahm4 following your desiderata) ! harmonic isopycnal or geopotential: ! ahm3 (ahm4) defined at u- (v-) point DO jk = 1, jpk DO jj = 2, jpj DO ji = 2, jpi ahm3(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji ,jj-1,jk) ) ahm4(ji,jj,jk) = 0.5 * ( ahm2(ji,jj,jk) + ahm2(ji-1,jj ,jk) ) END DO END DO END DO ahm3 ( :, 1, :) = ahm3 ( :, 2, :) ahm4 ( :, 1, :) = ahm4 ( :, 2, :) ! Lateral boundary conditions on ( ahm3, ahm4 ) ! ============== CALL lbc_lnk( ahm3, 'U', 1. ) ! U-point, unchanged sign CALL lbc_lnk( ahm4, 'V', 1. ) ! V-point, unchanged sign ! Control print IF( lwp .AND. ld_print ) THEN WRITE(numout,*) WRITE(numout,*) ' ahm3 array level 1' CALL prihre(ahm3(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) WRITE(numout,*) WRITE(numout,*) ' ahm4 array level 1' CALL prihre(ahm4(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) ENDIF END SUBROUTINE ldf_dyn_c3d_orca