Changeset 82 for trunk


Ignore:
Timestamp:
12/17/13 22:57:08 (10 years ago)
Author:
smasson
Message:

good viscosity and diffusivity in trop075

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

Legend:

Unmodified
Added
Removed
  • trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90

    r7 r82  
    2020   USE ldfslp          ! ??? 
    2121   USE ioipsl 
     22   USE iom 
    2223   USE in_out_manager  ! I/O manager 
    2324   USE lib_mpp         ! distribued memory computing library 
  • trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90

    r7 r82  
    7373         ! Special case for ORCA R1, R2 and R4 configurations (overwrite the value of ahm1 ahm2) 
    7474         ! ============================================== 
    75          IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) )   CALL ldf_dyn_c2d_orca( ld_print ) 
    76          IF( cp_cfg == "orca" .AND.   jp_cfg == 1)                       CALL ldf_dyn_c2d_orca_R1( ld_print ) 
     75         IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) )   CALL ldf_dyn_c2d_orca( ld_print, ld_cte = .true.  ) 
     76         IF( cp_cfg == "orca" .AND.   jp_cfg == 1                    )   CALL ldf_dyn_c2d_orca( ld_print, ld_cte = .false. ) 
     77         IF( cp_cfg == "trop" .AND.   jp_cfg == 075                  )   CALL ldf_dyn_c2d_orca( ld_print, ld_cte = .true.  ) 
    7778 
    7879         ! Control print 
     
    124125 
    125126 
    126    SUBROUTINE ldf_dyn_c2d_orca( ld_print ) 
     127   SUBROUTINE ldf_dyn_c2d_orca( ld_print, ld_cte ) 
    127128      !!---------------------------------------------------------------------- 
    128129      !!                 ***  ROUTINE ldf_dyn_c2d  *** 
     
    142143      ! 
    143144      LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout 
    144       ! 
    145       INTEGER  ::   ji, jj, jn   ! dummy loop indices 
    146       INTEGER  ::   inum, iim, ijm            ! local integers 
    147       INTEGER  ::   ifreq, il1, il2, ij, ii 
    148       REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk 
    149       CHARACTER (len=15) ::   clexp 
    150       INTEGER, POINTER, DIMENSION(:,:)  :: icof 
    151       INTEGER, POINTER, DIMENSION(:,:)  :: idata 
     145      LOGICAL, INTENT (in) ::   ld_cte     !  
     146      ! 
     147      INTEGER  ::   ji, jj   ! dummy loop indices 
     148      INTEGER  ::   inum            ! local integers 
     149      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk,zam20s 
     150      REAL(wp), POINTER, DIMENSION(:,:)  :: zcof 
    152151      !!---------------------------------------------------------------------- 
    153152      !                                 
    154       CALL wrk_alloc( jpi   , jpj   , icof  ) 
    155       CALL wrk_alloc( jpidta, jpjdta, idata ) 
     153      CALL wrk_alloc( jpi, jpj, zcof ) 
    156154      ! 
    157155      IF(lwp) WRITE(numout,*) 
     
    167165      ! Read 2d integer array to specify western boundary increase in the 
    168166      ! ===================== equatorial strip (20N-20S) defined at t-points 
    169  
    170       CALL ctl_opn( inum, 'ahmcoef', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
    171       READ(inum,9101) clexp, iim, ijm 
    172       READ(inum,'(/)') 
    173       ifreq = 40 
    174       il1 = 1 
    175       DO jn = 1, jpidta/ifreq+1 
    176          READ(inum,'(/)') 
    177          il2 = MIN( jpidta, il1+ifreq-1 ) 
    178          READ(inum,9201) ( ii, ji = il1, il2, 5 ) 
    179          READ(inum,'(/)') 
    180          DO jj = jpjdta, 1, -1 
    181             READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 ) 
    182          END DO 
    183          il1 = il1 + ifreq 
    184       END DO 
    185        
    186       DO jj = 1, nlcj 
    187          DO ji = 1, nlci 
    188             icof(ji,jj) = idata( mig(ji), mjg(jj) ) 
    189          END DO 
    190       END DO 
    191       DO jj = nlcj+1, jpj 
    192          DO ji = 1, nlci 
    193             icof(ji,jj) = icof(ji,nlcj) 
    194          END DO 
    195       END DO 
    196       DO jj = 1, jpj 
    197          DO ji = nlci+1, jpi 
    198             icof(ji,jj) = icof(nlci,jj) 
    199          END DO 
    200       END DO 
    201  
    202  9101 FORMAT(1x,a15,2i8) 
    203  9201 FORMAT(3x,13(i3,12x)) 
    204  9202 FORMAT(i3,41i3) 
    205  
     167       
     168      CALL iom_open ( 'ahmcoef.nc', inum )   
     169      CALL iom_get  ( inum, jpdom_data, 'ahmcoef', zcof ) 
     170      CALL iom_close( inum ) 
    206171 
    207172      ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator) 
     
    218183 
    219184      zahmeq = aht0 
    220        
     185      IF (ld_cte) THEN   ;   zam20s = ahm0 
     186      ELSE               ;   zam20s = ahm0 * COS( rad * 20. ) 
     187      END IF 
     188     
    221189      DO jj = 1, jpj 
    222190         DO ji = 1, jpi 
    223191            IF( ABS( gphif(ji,jj) ) >= 20. ) THEN 
    224                ahm2(ji,jj) = ahm0 
     192               IF (ld_cte) ahm2(ji,jj) = ahm0 
    225193            ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN 
    226                ahm2(ji,jj) =  zahmeq 
     194               ahm2(ji,jj) = zahmeq 
    227195            ELSE 
    228                ahm2(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   & 
     196               ahm2(ji,jj) = zahmeq + (zam20s-zahmeq)/2.   & 
    229197                  * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) ) 
    230198            ENDIF 
    231199            IF( ABS( gphit(ji,jj) ) >= 20. ) THEN 
    232                ahm1(ji,jj) = ahm0 
     200               IF (ld_cte) ahm1(ji,jj) = ahm0 
    233201            ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN 
    234                ahm1(ji,jj) =  zahmeq 
     202               ahm1(ji,jj) = zahmeq 
    235203            ELSE 
    236                ahm1(ji,jj) = zahmeq + (ahm0-zahmeq)/2.   & 
     204               ahm1(ji,jj) = zahmeq + (zam20s-zahmeq)/2.   & 
    237205                  * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) ) 
    238206            ENDIF 
     
    244212      DO jj = 1, jpjm1 
    245213         DO ji = 1, jpim1 
    246             zcoft = FLOAT( icof(ji,jj) ) / 100. 
     214            zcoft = zcof(ji,jj) / 100. 
    247215            ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)  
    248216         END DO 
    249217      END DO 
    250218      ! f-point 
    251       icof(:,:) = icof(:,:) * tmask(:,:,1) 
     219      zcof(:,:) = zcof(:,:) * tmask(:,:,1) 
    252220      DO jj = 1, jpjm1 
    253221         DO ji = 1, jpim1   ! NO vector opt. 
     
    256224               zcoff = 1. 
    257225            ELSE 
    258                zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   & 
     226               zcoff = ( zcof(ji,jj+1) + zcof(ji+1,jj+1) + zcof(ji,jj) + zcof(ji,jj+1) )   & 
    259227                     / (zmsk * 100.) 
    260228            ENDIF 
     
    279247      ENDIF 
    280248      ! 
    281       CALL wrk_dealloc( jpi   , jpj   , icof  ) 
    282       CALL wrk_dealloc( jpidta, jpjdta, idata ) 
     249      CALL wrk_dealloc( jpi, jpj, zcof ) 
    283250      ! 
    284251   END SUBROUTINE ldf_dyn_c2d_orca 
    285252 
    286  
    287    SUBROUTINE ldf_dyn_c2d_orca_R1( ld_print ) 
    288       !!---------------------------------------------------------------------- 
    289       !!                 ***  ROUTINE ldf_dyn_c2d  *** 
    290       !! 
    291       !!                   **** W A R N I N G **** 
    292       !! 
    293       !!                ORCA R1 configuration 
    294       !!                   
    295       !!                   **** W A R N I N G **** 
    296       !!                   
    297       !! ** Purpose :   initializations of the lateral viscosity for orca R1 
    298       !! 
    299       !! ** Method  :   blah blah blah... 
    300       !! 
    301       !!---------------------------------------------------------------------- 
    302       USE ldftra_oce, ONLY:   aht0 
    303       ! 
    304       LOGICAL, INTENT (in) ::   ld_print   ! If true, output arrays on numout 
    305       ! 
    306       INTEGER ::   ji, jj, jn      ! dummy loop indices 
    307       INTEGER ::   inum            ! temporary logical unit 
    308       INTEGER ::   iim, ijm 
    309       INTEGER ::   ifreq, il1, il2, ij, ii 
    310       REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk, zam20s 
    311       CHARACTER (len=15) ::   clexp 
    312       INTEGER, POINTER, DIMENSION(:,:)  :: icof 
    313       INTEGER, POINTER, DIMENSION(:,:)  :: idata 
    314       !!---------------------------------------------------------------------- 
    315       !                                 
    316       CALL wrk_alloc( jpi   , jpj   , icof  ) 
    317       CALL wrk_alloc( jpidta, jpjdta, idata ) 
    318       !                                 
    319  
    320       IF(lwp) WRITE(numout,*) 
    321       IF(lwp) WRITE(numout,*) 'inildf: 2d eddy viscosity coefficient' 
    322       IF(lwp) WRITE(numout,*) '~~~~~~  --' 
    323       IF(lwp) WRITE(numout,*) '        orca_r1 configuration' 
    324  
    325 #if defined key_antarctic 
    326 #     include "ldfdyn_antarctic.h90" 
    327 #elif defined key_arctic 
    328 #     include "ldfdyn_arctic.h90" 
    329 #else 
    330       ! Read 2d integer array to specify western boundary increase in the 
    331       ! ===================== equatorial strip (20N-20S) defined at t-points 
    332  
    333       CALL ctl_opn( inum, 'ahmcoef', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   & 
    334          &           1, numout, lwp ) 
    335       REWIND inum 
    336       READ(inum,9101) clexp, iim, ijm 
    337       READ(inum,'(/)') 
    338       ifreq = 40 
    339       il1 = 1 
    340       DO jn = 1, jpidta/ifreq+1 
    341          READ(inum,'(/)') 
    342          il2 = MIN( jpidta, il1+ifreq-1 ) 
    343          READ(inum,9201) ( ii, ji = il1, il2, 5 ) 
    344          READ(inum,'(/)') 
    345          DO jj = jpjdta, 1, -1 
    346             READ(inum,9202) ij, ( idata(ji,jj), ji = il1, il2 ) 
    347          END DO 
    348          il1 = il1 + ifreq 
    349       END DO 
    350        
    351       DO jj = 1, nlcj 
    352          DO ji = 1, nlci 
    353             icof(ji,jj) = idata( mig(ji), mjg(jj) ) 
    354          END DO 
    355       END DO 
    356       DO jj = nlcj+1, jpj 
    357          DO ji = 1, nlci 
    358             icof(ji,jj) = icof(ji,nlcj) 
    359          END DO 
    360       END DO 
    361       DO jj = 1, jpj 
    362          DO ji = nlci+1, jpi 
    363             icof(ji,jj) = icof(nlci,jj) 
    364          END DO 
    365       END DO 
    366  
    367  9101 FORMAT(1x,a15,2i8) 
    368  9201 FORMAT(3x,13(i3,12x)) 
    369  9202 FORMAT(i3,41i3) 
    370  
    371  
    372       ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator) 
    373       ! ================= 
    374       ! define ahm1 and ahm2 at the right grid point position 
    375       ! (USER: modify ahm1 and ahm2 following your desiderata) 
    376        
    377        
    378       ! Decrease ahm to zahmeq m2/s in the tropics 
    379       ! (from 90   to 20   degrees: ahm = scaled by local metrics 
    380       !  from 20   to  2.5 degrees: ahm = decrease in (1-cos)/2 
    381       !  from  2.5 to  0   degrees: ahm = constant 
    382       ! symmetric in the south hemisphere) 
    383  
    384       zahmeq = aht0 
    385       zam20s = ahm0*COS( rad * 20. ) 
    386        
    387       DO jj = 1, jpj 
    388          DO ji = 1, jpi 
    389             IF( ABS( gphif(ji,jj) ) >= 20. ) THEN 
    390 !              leave as set in ldf_dyn_c2d 
    391             ELSEIF( ABS( gphif(ji,jj) ) <= 2.5 ) THEN 
    392                ahm2(ji,jj) =  zahmeq 
    393             ELSE 
    394                ahm2(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   & 
    395                   * ( 1. - COS( rad * ( ABS(gphif(ji,jj))-2.5 ) * 180. / 17.5 ) ) 
    396             ENDIF 
    397             IF( ABS( gphit(ji,jj) ) >= 20. ) THEN 
    398 !             leave as set in ldf_dyn_c2d 
    399             ELSEIF( ABS( gphit(ji,jj) ) <= 2.5 ) THEN 
    400                ahm1(ji,jj) =  zahmeq 
    401             ELSE 
    402                ahm1(ji,jj) =  zahmeq + (zam20s-zahmeq)/2.   & 
    403                   * ( 1. - COS( rad * ( ABS(gphit(ji,jj))-2.5 ) * 180. / 17.5 ) ) 
    404             ENDIF 
    405          END DO 
    406       END DO 
    407  
    408       ! increase along western boundaries of equatorial strip 
    409       ! t-point 
    410       DO jj = 1, jpjm1 
    411          DO ji = 1, jpim1 
    412           IF( ABS( gphit(ji,jj) ) < 20. ) THEN 
    413             zcoft = FLOAT( icof(ji,jj) ) / 100. 
    414             ahm1(ji,jj) = zcoft * ahm0 + (1.-zcoft) * ahm1(ji,jj)  
    415           ENDIF 
    416          END DO 
    417       END DO 
    418       ! f-point 
    419       icof(:,:) = icof(:,:) * tmask(:,:,1) 
    420       DO jj = 1, jpjm1 
    421          DO ji = 1, jpim1 
    422           IF( ABS( gphif(ji,jj) ) < 20. ) THEN 
    423             zmsk = tmask(ji,jj+1,1) + tmask(ji+1,jj+1,1) + tmask(ji,jj,1) + tmask(ji,jj+1,1) 
    424             IF( zmsk == 0. ) THEN 
    425                zcoff = 1. 
    426             ELSE 
    427                zcoff = FLOAT( icof(ji,jj+1) + icof(ji+1,jj+1) + icof(ji,jj) + icof(ji,jj+1) )   & 
    428                      / (zmsk * 100.) 
    429             ENDIF 
    430             ahm2(ji,jj) = zcoff * ahm0 + (1.-zcoff) * ahm2(ji,jj) 
    431           ENDIF 
    432          END DO 
    433       END DO 
    434 #endif 
    435        
    436       ! Lateral boundary conditions on ( ahm1, ahm2 ) 
    437       !                                ============== 
    438       CALL lbc_lnk( ahm1, 'T', 1. )   ! T-point, unchanged sign 
    439       CALL lbc_lnk( ahm2, 'F', 1. )   ! F-point, unchanged sign 
    440  
    441       ! Control print 
    442       IF( lwp .AND. ld_print ) THEN 
    443          WRITE(numout,*) 
    444          WRITE(numout,*) 'inildf: 2D ahm1 array' 
    445          CALL prihre(ahm1,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    446          WRITE(numout,*) 
    447          WRITE(numout,*) 'inildf: 2D ahm2 array' 
    448          CALL prihre(ahm2,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    449       ENDIF 
    450       ! 
    451       CALL wrk_dealloc( jpi   , jpj   , icof  ) 
    452       CALL wrk_dealloc( jpidta, jpjdta, idata ) 
    453       ! 
    454    END SUBROUTINE ldf_dyn_c2d_orca_R1 
  • trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90

    r7 r82  
    188188 
    189189      ! ORCA R05: Take the minimum between aeiw  and aeiv0 
    190       IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN 
     190      IF( ( cp_cfg == "orca" .AND. jp_cfg == 05 ) .OR. ( cp_cfg == "trop" .AND. jp_cfg == 075 ) ) THEN 
    191191         DO jj = 2, jpjm1 
    192192            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    214214 
    215215      ! ORCA R05: add a space variation on aht (=aeiv except at the equator and river mouth) 
    216       IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN 
     216      IF( ( cp_cfg == "orca" .AND. jp_cfg == 05 ) .OR. ( cp_cfg == "trop" .AND. jp_cfg == 075 ) ) THEN 
    217217         zf20     = 2._wp * omega * SIN( rad * 20._wp ) 
    218218         zaht_min = 100._wp                           ! minimum value for aht 
  • trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_c2d.h90

    r1 r82  
    7171         ! Special case for ORCA R2 and R4 configurations (overwrite the value of ahtu ahtv ahtw) 
    7272         ! ============================================== 
    73          IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) )   THEN 
     73         IF(     ( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) ) & 
     74            .OR. ( cp_cfg == "trop" .AND.   jp_cfg == 075                  ) )  THEN 
    7475            ahtu(:,:) = aht0              ! set ahtu = ahtv at u- and v-points, 
    7576            ahtv(:,:) = aht0              ! and ahtw at w-point 
Note: See TracChangeset for help on using the changeset viewer.