!!---------------------------------------------------------------------- !! *** ldfdyn_c2d.h90 *** !!---------------------------------------------------------------------- !! ldf_dyn_c2d : set the lateral viscosity coefficients !! ldf_dyn_c2d_orca : specific case for orca r2 and r4 !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! OPA 9.0 , LOCEAN-IPSL (2005) !! $Id$ !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt !!---------------------------------------------------------------------- SUBROUTINE ldf_dyn_c2d( ld_print ) !!---------------------------------------------------------------------- !! *** ROUTINE ldf_dyn_c2d *** !! !! ** Purpose : initializations of the horizontal ocean physics !! !! ** Method : !! 2D eddy viscosity coefficients ( longitude, latitude ) !! !! harmonic operator : ahm1 is defined at t-point !! ahm2 is defined at f-point !! + isopycnal : ahm3 is defined at u-point !! or geopotential ahm4 is defined at v-point !! iso-model level : ahm3, ahm4 not used !! !! biharmonic operator : ahm3 is defined at u-point !! ahm4 is defined at v-point !! : ahm1, ahm2 not used !! !!---------------------------------------------------------------------- !! * Arguments LOGICAL, INTENT (in) :: ld_print ! If true, output arrays on numout !! * Local variables INTEGER :: ji, jj REAL(wp) :: za00, zd_max, zetmax, zeumax, zefmax, zevmax !!---------------------------------------------------------------------- IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'ldf_dyn_c2d : 2d lateral eddy viscosity coefficient' IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' IF(lwp) WRITE(numout,*) ! harmonic operator (ahm1, 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 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) = za00 * zetmax ahm2(ji,jj) = za00 * zefmax END DO END DO 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 ! 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 ) ) CALL ldf_dyn_c2d_orca( ld_print ) ! Control print IF( lwp .AND. ld_print ) THEN WRITE(numout,*) WRITE(numout,*) 'inildf: 2D ahm1 array' CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) WRITE(numout,*) WRITE(numout,*) 'inildf: 2D ahm2 array' CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) ENDIF ENDIF ! biharmonic operator (ahm3, ahm4) : at U- and V-points (used for bilaplacian operator ! ================================= whatever its orientation is) IF( ln_dynldf_bilap ) THEN ! (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 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) = za00 * zeumax * zeumax * zeumax ahm4(ji,jj) = za00 * zevmax * zevmax * zevmax END DO END DO ! Control print IF( lwp .AND. ld_print ) THEN WRITE(numout,*) WRITE(numout,*) 'inildf: ahm3 array' CALL prihre(ahm3,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) WRITE(numout,*) WRITE(numout,*) 'inildf: ahm4 array' CALL prihre(ahm4,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) ENDIF ENDIF END SUBROUTINE ldf_dyn_c2d SUBROUTINE ldf_dyn_c2d_orca( ld_print ) !!---------------------------------------------------------------------- !! *** ROUTINE ldf_dyn_c2d *** !! !! **** W A R N I N G **** !! !! ORCA R2 and R4 configurations !! !! **** W A R N I N G **** !! !! ** Purpose : initializations of the lateral viscosity for orca R2 !! !! ** 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, jn ! dummy loop indices INTEGER :: inum ! temporary logical unit INTEGER :: iim, ijm INTEGER :: ifreq, il1, il2, ij, ii INTEGER, DIMENSION(jpidta,jpidta) :: idata INTEGER, DIMENSION(jpi ,jpj ) :: icof REAL(wp) :: zahmeq, zcoft, zcoff, zmsk CHARACTER (len=15) :: clexp !!---------------------------------------------------------------------- IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient' IF(lwp) WRITE(numout,*) '~~~~~~ --' IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' orca ocean model' IF(lwp) WRITE(numout,*) #if defined key_antarctic # include "ldfdyn_antarctic.h90" #elif defined key_arctic # include "ldfdyn_arctic.h90" #else ! 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 ( T- and F- points) (used for laplacian operator) ! ================= ! define ahm1 and ahm2 at the right grid point position ! (USER: modify ahm1 and ahm2 following your desiderata) ! 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) zahmeq = aht0 DO jj = 1, jpj DO ji = 1, jpi IF( ABS( gphif(ji,jj) ) >= 20. ) THEN ahm2(ji,jj) = ahm0 ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN ahm2(ji,jj) = zahmeq ELSE ahm2(ji,jj) = zahmeq + (ahm0-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) = ahm0 ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN ahm1(ji,jj) = zahmeq ELSE ahm1(ji,jj) = zahmeq + (ahm0-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) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj) 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) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj) END DO 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 .AND. ld_print ) THEN WRITE(numout,*) WRITE(numout,*) 'inildf: 2D ahm1 array' CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) WRITE(numout,*) WRITE(numout,*) 'inildf: 2D ahm2 array' CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) ENDIF END SUBROUTINE ldf_dyn_c2d_orca