Ignore:
Timestamp:
2015-07-15T17:46:12+02:00 (5 years ago)
Author:
andrewryan
Message:

merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90

    r5034 r5600  
    2121   !!   tra_dmp       : update the tracer trend with the internal damping 
    2222   !!   tra_dmp_init  : initialization, namlist read, parameters control 
    23    !!   dtacof_zoom   : restoring coefficient for zoom domain 
    24    !!   dtacof        : restoring coefficient for global domain 
    25    !!   cofdis        : compute the distance to the coastline 
    2623   !!---------------------------------------------------------------------- 
    2724   USE oce            ! ocean: variables 
     
    3936   USE wrk_nemo       ! Memory allocation 
    4037   USE timing         ! Timing 
     38   USE iom 
    4139 
    4240   IMPLICIT NONE 
     
    4543   PUBLIC   tra_dmp      ! routine called by step.F90 
    4644   PUBLIC   tra_dmp_init ! routine called by opa.F90 
    47    PUBLIC   dtacof       ! routine called by tradmp.F90, trcdmp.F90 and dyndmp.F90 
    48    PUBLIC   dtacof_zoom  ! routine called by tradmp.F90, trcdmp.F90 and dyndmp.F90 
    49  
    50 !!gm  why all namelist variable public????   only ln_tradmp should be sufficient 
    5145 
    5246   !                               !!* Namelist namtra_dmp : T & S newtonian damping * 
     47   ! nn_zdmp and cn_resto are public as they are used by C1D/dyndmp.F90 
    5348   LOGICAL , PUBLIC ::   ln_tradmp   !: internal damping flag 
    54    INTEGER , PUBLIC ::   nn_hdmp     ! = 0/-1/'latitude' for damping over T and S 
    5549   INTEGER , PUBLIC ::   nn_zdmp     ! = 0/1/2 flag for damping in the mixed layer 
    56    REAL(wp), PUBLIC ::   rn_surf     ! surface time scale for internal damping        [days] 
    57    REAL(wp), PUBLIC ::   rn_bot      ! bottom time scale for internal damping         [days] 
    58    REAL(wp), PUBLIC ::   rn_dep      ! depth of transition between rn_surf and rn_bot [meters] 
    59    INTEGER , PUBLIC ::   nn_file     ! = 1 create a damping.coeff NetCDF file  
     50   CHARACTER(LEN=200) , PUBLIC :: cn_resto      ! name of netcdf file containing restoration coefficient field 
     51   ! 
     52 
    6053 
    6154   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   strdmp   !: damping salinity trend (psu/s) 
     
    197190      !! ** Method  :   read the namtra_dmp namelist and check the parameters 
    198191      !!---------------------------------------------------------------------- 
    199       INTEGER  ::   ios   ! Local integer output status for namelist read 
    200       !! 
    201       NAMELIST/namtra_dmp/ ln_tradmp, nn_hdmp, nn_zdmp, rn_surf, rn_bot, rn_dep, nn_file 
    202       !!---------------------------------------------------------------------- 
    203       ! 
    204       REWIND( numnam_ref )              ! Namelist namtra_dmp in reference namelist : Temperature and salinity damping term 
     192      NAMELIST/namtra_dmp/ ln_tradmp, nn_zdmp, cn_resto 
     193      INTEGER ::  ios         ! Local integer for output status of namelist read 
     194      INTEGER :: imask        ! File handle  
     195      !! 
     196      !!---------------------------------------------------------------------- 
     197      ! 
     198      REWIND( numnam_ref )   ! Namelist namtra_dmp in reference namelist : T & S relaxation 
    205199      READ  ( numnam_ref, namtra_dmp, IOSTAT = ios, ERR = 901) 
    206200901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in reference namelist', lwp ) 
    207201      ! 
    208       REWIND( numnam_cfg )              ! Namelist namtra_dmp in configuration namelist : Temperature and salinity damping term 
     202      REWIND( numnam_cfg )   ! Namelist namtra_dmp in configuration namelist : T & S relaxation 
    209203      READ  ( numnam_cfg, namtra_dmp, IOSTAT = ios, ERR = 902 ) 
    210204902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_dmp in configuration namelist', lwp ) 
    211205      IF(lwm) WRITE ( numond, namtra_dmp ) 
    212        
    213       IF( lzoom .AND. .NOT. lk_c1d )   nn_zdmp = 0          ! restoring to climatology at closed north or south boundaries 
    214  
    215       IF(lwp) THEN                       ! Namelist print 
     206 
     207      IF(lwp) THEN                 !Namelist print 
    216208         WRITE(numout,*) 
    217          WRITE(numout,*) 'tra_dmp_init : T and S newtonian damping' 
     209         WRITE(numout,*) 'tra_dmp_init : T and S newtonian relaxation' 
    218210         WRITE(numout,*) '~~~~~~~' 
    219          WRITE(numout,*) '   Namelist namtra_dmp : set damping parameter' 
    220          WRITE(numout,*) '      add a damping term or not       ln_tradmp = ', ln_tradmp 
    221          WRITE(numout,*) '      T and S damping option          nn_hdmp   = ', nn_hdmp 
    222          WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp, '(non-C1D zoom: forced to 0)' 
    223          WRITE(numout,*) '      surface time scale (days)       rn_surf   = ', rn_surf 
    224          WRITE(numout,*) '      bottom time scale (days)        rn_bot    = ', rn_bot 
    225          WRITE(numout,*) '      depth of transition (meters)    rn_dep    = ', rn_dep 
    226          WRITE(numout,*) '      create a damping.coeff file     nn_file   = ', nn_file 
     211         WRITE(numout,*) '   Namelist namtra_dmp : set relaxation parameters' 
     212         WRITE(numout,*) '      Apply relaxation   or not       ln_tradmp = ', ln_tradmp 
     213         WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp 
     214         WRITE(numout,*) '      Damping file name               cn_resto  = ', cn_resto 
    227215         WRITE(numout,*) 
    228216      ENDIF 
    229217 
    230       IF( ln_tradmp ) THEN               ! initialization for T-S damping 
    231          ! 
     218      IF( ln_tradmp) THEN 
     219         ! 
     220         !Allocate arrays 
    232221         IF( tra_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_dmp_init: unable to allocate arrays' ) 
    233          ! 
    234 !!gm  I don't understand the specificities of c1d case......    
    235 !!gm  to be check with the autor of these lines 
    236           
    237 #if ! defined key_c1d 
    238          SELECT CASE ( nn_hdmp ) 
    239          CASE (  -1  )   ;   IF(lwp) WRITE(numout,*) '   tracer damping in the Med & Red seas only' 
    240          CASE ( 1:90 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping poleward of', nn_hdmp, ' degrees' 
    241          CASE DEFAULT 
    242             WRITE(ctmp1,*) '          bad flag value for nn_hdmp = ', nn_hdmp 
    243             CALL ctl_stop(ctmp1) 
     222 
     223         !Check values of nn_zdmp 
     224         SELECT CASE (nn_zdmp) 
     225         CASE ( 0 )  ; IF(lwp) WRITE(numout,*) '   tracer damping as specified by mask' 
     226         CASE ( 1 )  ; IF(lwp) WRITE(numout,*) '   no tracer damping in the turbocline' 
     227         CASE ( 2 )  ; IF(lwp) WRITE(numout,*) '   no tracer damping in the mixed layer' 
    244228         END SELECT 
    245          ! 
    246 #endif 
    247          SELECT CASE ( nn_zdmp ) 
    248          CASE ( 0 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping throughout the water column' 
    249          CASE ( 1 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the turbocline (avt > 5 cm2/s)' 
    250          CASE ( 2 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the mixed layer' 
    251          CASE DEFAULT 
    252             WRITE(ctmp1,*) 'bad flag value for nn_zdmp = ', nn_zdmp 
    253             CALL ctl_stop(ctmp1) 
    254          END SELECT 
    255          ! 
     229 
     230         !TG: Initialisation of dtatsd - Would it be better to have dmpdta routine 
     231         !so can damp to something other than intitial conditions files? 
    256232         IF( .NOT.ln_tsd_tradmp ) THEN 
    257233            CALL ctl_warn( 'tra_dmp_init: read T-S data not initialized, we force ln_tsd_tradmp=T' ) 
    258234            CALL dta_tsd_init( ld_tradmp=ln_tradmp )        ! forces the initialisation of T-S data 
    259235         ENDIF 
    260          ! 
    261          strdmp(:,:,:) = 0._wp       ! internal damping salinity trend (used in asmtrj) 
     236 
     237         !initialise arrays - Are these actually used anywhere else? 
     238         strdmp(:,:,:) = 0._wp 
    262239         ttrdmp(:,:,:) = 0._wp 
    263          !                          ! Damping coefficients initialization 
    264          IF( lzoom .AND. .NOT. lk_c1d ) THEN   ;   CALL dtacof_zoom( resto ) 
    265          ELSE               ;   CALL dtacof( nn_hdmp, rn_surf, rn_bot, rn_dep, nn_file, 'TRA', resto ) 
    266          ENDIF 
    267          ! 
    268       ENDIF 
    269       ! 
     240 
     241         !Read in mask from file 
     242         CALL iom_open ( cn_resto, imask) 
     243         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto) 
     244         CALL iom_close( imask ) 
     245       ENDIF 
     246 
    270247   END SUBROUTINE tra_dmp_init 
    271248 
    272  
    273    SUBROUTINE dtacof_zoom( presto ) 
    274       !!---------------------------------------------------------------------- 
    275       !!                  ***  ROUTINE dtacof_zoom  *** 
    276       !! 
    277       !! ** Purpose :   Compute the damping coefficient for zoom domain 
    278       !! 
    279       !! ** Method  : - set along closed boundary due to zoom a damping over 
    280       !!                6 points with a max time scale of 5 days. 
    281       !!              - ORCA arctic/antarctic zoom: set the damping along 
    282       !!                south/north boundary over a latitude strip. 
    283       !! 
    284       !! ** Action  : - resto, the damping coeff. for T and S 
    285       !!---------------------------------------------------------------------- 
    286       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)  ::   presto   ! restoring coeff. (s-1) 
    287       ! 
    288       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    289       REAL(wp) ::   zlat, zlat0, zlat1, zlat2, z1_5d   ! local scalar 
    290       REAL(wp), DIMENSION(6)  ::   zfact               ! 1Dworkspace 
    291       !!---------------------------------------------------------------------- 
    292       ! 
    293       IF( nn_timing == 1 )  CALL timing_start( 'dtacof_zoom') 
    294       ! 
    295  
    296       zfact(1) =  1._wp 
    297       zfact(2) =  1._wp 
    298       zfact(3) = 11._wp / 12._wp 
    299       zfact(4) =  8._wp / 12._wp 
    300       zfact(5) =  4._wp / 12._wp 
    301       zfact(6) =  1._wp / 12._wp 
    302       zfact(:) = zfact(:) / ( 5._wp * rday )    ! 5 days max restoring time scale 
    303  
    304       presto(:,:,:) = 0._wp 
    305  
    306       ! damping along the forced closed boundary over 6 grid-points 
    307       DO jn = 1, 6 
    308          IF( lzoom_w )   presto( mi0(jn+jpizoom):mi1(jn+jpizoom), : , : )                    = zfact(jn)   ! west  closed 
    309          IF( lzoom_s )   presto( : , mj0(jn+jpjzoom):mj1(jn+jpjzoom), : )                    = zfact(jn)   ! south closed  
    310          IF( lzoom_e )   presto( mi0(jpiglo+jpizoom-1-jn):mi1(jpiglo+jpizoom-1-jn) , : , : ) = zfact(jn)   ! east  closed  
    311          IF( lzoom_n )   presto( : , mj0(jpjglo+jpjzoom-1-jn):mj1(jpjglo+jpjzoom-1-jn) , : ) = zfact(jn)   ! north closed 
    312       END DO 
    313  
    314       !                                           ! ==================================================== 
    315       IF( cp_cfz == "arctic" .OR. cp_cfz == "antarctic" ) THEN   !  ORCA configuration : arctic or antarctic zoom 
    316          !                                        ! ==================================================== 
    317          IF(lwp) WRITE(numout,*) 
    318          IF(lwp .AND. cp_cfz == "arctic" ) WRITE(numout,*) '              dtacof_zoom : ORCA    Arctic zoom' 
    319          IF(lwp .AND. cp_cfz == "antarctic" ) WRITE(numout,*) '           dtacof_zoom : ORCA Antarctic zoom' 
    320          IF(lwp) WRITE(numout,*) 
    321          ! 
    322          !                          ! Initialization :  
    323          presto(:,:,:) = 0._wp 
    324          zlat0 = 10._wp                     ! zlat0 : latitude strip where resto decreases 
    325          zlat1 = 30._wp                     ! zlat1 : resto = 1 before zlat1 
    326          zlat2 = zlat1 + zlat0              ! zlat2 : resto decreases from 1 to 0 between zlat1 and zlat2 
    327          z1_5d = 1._wp / ( 5._wp * rday )   ! z1_5d : 1 / 5days 
    328  
    329          DO jk = 2, jpkm1           ! Compute arrays resto ; value for internal damping : 5 days 
    330             DO jj = 1, jpj 
    331                DO ji = 1, jpi 
    332                   zlat = ABS( gphit(ji,jj) ) 
    333                   IF( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN 
    334                      presto(ji,jj,jk) = 0.5_wp * z1_5d * (  1._wp - COS( rpi*(zlat2-zlat)/zlat0 )  )  
    335                   ELSEIF( zlat < zlat1 ) THEN 
    336                      presto(ji,jj,jk) = z1_5d 
    337                   ENDIF 
    338                END DO 
    339             END DO 
    340          END DO 
    341          ! 
    342       ENDIF 
    343       !                             ! Mask resto array 
    344       presto(:,:,:) = presto(:,:,:) * tmask(:,:,:) 
    345       ! 
    346       IF( nn_timing == 1 )  CALL timing_stop( 'dtacof_zoom') 
    347       ! 
    348    END SUBROUTINE dtacof_zoom 
    349  
    350  
    351    SUBROUTINE dtacof( kn_hdmp, pn_surf, pn_bot, pn_dep,  & 
    352       &               kn_file, cdtype , presto           ) 
    353       !!---------------------------------------------------------------------- 
    354       !!                  ***  ROUTINE dtacof  *** 
    355       !! 
    356       !! ** Purpose :   Compute the damping coefficient 
    357       !! 
    358       !! ** Method  :   Arrays defining the damping are computed for each grid 
    359       !!                point for temperature and salinity (resto) 
    360       !!                Damping depends on distance to coast, depth and latitude 
    361       !! 
    362       !! ** Action  : - resto, the damping coeff. for T and S 
    363       !!---------------------------------------------------------------------- 
    364       USE iom 
    365       USE ioipsl 
    366       !! 
    367       INTEGER                         , INTENT(in   )  ::  kn_hdmp    ! damping option 
    368       REAL(wp)                        , INTENT(in   )  ::  pn_surf    ! surface time scale (days) 
    369       REAL(wp)                        , INTENT(in   )  ::  pn_bot     ! bottom time scale (days) 
    370       REAL(wp)                        , INTENT(in   )  ::  pn_dep     ! depth of transition (meters) 
    371       INTEGER                         , INTENT(in   )  ::  kn_file    ! save the damping coef on a file or not 
    372       CHARACTER(len=3)                , INTENT(in   )  ::  cdtype     ! =TRA, TRC or DYN (tracer/dynamics indicator) 
    373       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)  ::  presto     ! restoring coeff. (s-1) 
    374       ! 
    375       INTEGER  ::   ji, jj, jk                  ! dummy loop indices 
    376       INTEGER  ::   ii0, ii1, ij0, ij1          ! local integers 
    377       INTEGER  ::   inum0, icot                 !   -       - 
    378       REAL(wp) ::   zinfl, zlon                 ! local scalars 
    379       REAL(wp) ::   zlat, zlat0, zlat1, zlat2   !   -      - 
    380       REAL(wp) ::   zsdmp, zbdmp                !   -      - 
    381       CHARACTER(len=20)                   :: cfile 
    382       REAL(wp), POINTER, DIMENSION(:    ) :: zhfac  
    383       REAL(wp), POINTER, DIMENSION(:,:  ) :: zmrs  
    384       REAL(wp), POINTER, DIMENSION(:,:,:) :: zdct  
    385       !!---------------------------------------------------------------------- 
    386       ! 
    387       IF( nn_timing == 1 )  CALL timing_start('dtacof') 
    388       ! 
    389       CALL wrk_alloc( jpk, zhfac          ) 
    390       CALL wrk_alloc( jpi, jpj, zmrs      ) 
    391       CALL wrk_alloc( jpi, jpj, jpk, zdct ) 
    392 #if defined key_c1d 
    393       !                                   ! ==================== 
    394       !                                   !  C1D configuration : local domain 
    395       !                                   ! ==================== 
    396       ! 
    397       IF(lwp) WRITE(numout,*) 
    398       IF(lwp) WRITE(numout,*) '              dtacof : C1D 3x3 local domain' 
    399       IF(lwp) WRITE(numout,*) '              -----------------------------' 
    400       ! 
    401       presto(:,:,:) = 0._wp 
    402       ! 
    403       zsdmp = 1._wp / ( pn_surf * rday ) 
    404       zbdmp = 1._wp / ( pn_bot  * rday ) 
    405       DO jk = 2, jpkm1 
    406          DO jj = 1, jpj 
    407             DO ji = 1, jpi 
    408                !   ONLY vertical variation from zsdmp (sea surface) to zbdmp (bottom) 
    409                presto(ji,jj,jk) = zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(ji,jj,jk)/pn_dep) 
    410             END DO 
    411          END DO 
    412       END DO 
    413       ! 
    414       presto(:,:, : ) = presto(:,:,:) * tmask(:,:,:) 
    415 #else 
    416       !                                   ! ==================== 
    417       !                                   !  ORCA configuration : global domain 
    418       !                                   ! ==================== 
    419       ! 
    420       IF(lwp) WRITE(numout,*) 
    421       IF(lwp) WRITE(numout,*) '              dtacof : Global domain of ORCA' 
    422       IF(lwp) WRITE(numout,*) '              ------------------------------' 
    423       ! 
    424       presto(:,:,:) = 0._wp 
    425       ! 
    426       IF( kn_hdmp > 0 ) THEN      !  Damping poleward of 'nn_hdmp' degrees  ! 
    427          !                        !-----------------------------------------! 
    428          IF(lwp) WRITE(numout,*) 
    429          IF(lwp) WRITE(numout,*) '              Damping poleward of ', kn_hdmp, ' deg.' 
    430          ! 
    431          CALL iom_open ( 'dist.coast.nc', icot, ldstop = .FALSE. ) 
    432          ! 
    433          IF( icot > 0 ) THEN          ! distance-to-coast read in file 
    434             CALL iom_get  ( icot, jpdom_data, 'Tcoast', zdct ) 
    435             CALL iom_close( icot ) 
    436          ELSE                         ! distance-to-coast computed and saved in file (output in zdct) 
    437             CALL cofdis( zdct ) 
    438          ENDIF 
    439  
    440          !                            ! Compute arrays resto  
    441          zinfl = 1000.e3_wp                ! distance of influence for damping term 
    442          zlat0 = 10._wp                    ! latitude strip where resto decreases 
    443          zlat1 = REAL( kn_hdmp )           ! resto = 0 between -zlat1 and zlat1 
    444          zlat2 = zlat1 + zlat0             ! resto increases from 0 to 1 between |zlat1| and |zlat2| 
    445  
    446          DO jj = 1, jpj 
    447             DO ji = 1, jpi 
    448                zlat = ABS( gphit(ji,jj) ) 
    449                IF ( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN 
    450                   presto(ji,jj,1) = 0.5_wp * (  1._wp - COS( rpi*(zlat-zlat1)/zlat0 )  ) 
    451                ELSEIF ( zlat > zlat2 ) THEN 
    452                   presto(ji,jj,1) = 1._wp 
    453                ENDIF 
    454             END DO 
    455          END DO 
    456  
    457          IF ( kn_hdmp == 20 ) THEN       ! North Indian ocean (20N/30N x 45E/100E) : resto=0 
    458             DO jj = 1, jpj 
    459                DO ji = 1, jpi 
    460                   zlat = gphit(ji,jj) 
    461                   zlon = MOD( glamt(ji,jj), 360._wp ) 
    462                   IF ( zlat1 < zlat .AND. zlat < zlat2 .AND. 45._wp < zlon .AND. zlon < 100._wp ) THEN 
    463                      presto(ji,jj,1) = 0._wp 
    464                   ENDIF 
    465                END DO 
    466             END DO 
    467          ENDIF 
    468  
    469          zsdmp = 1._wp / ( pn_surf * rday ) 
    470          zbdmp = 1._wp / ( pn_bot  * rday ) 
    471          DO jk = 2, jpkm1 
    472             DO jj = 1, jpj 
    473                DO ji = 1, jpi 
    474                   zdct(ji,jj,jk) = MIN( zinfl, zdct(ji,jj,jk) ) 
    475                   !   ... Decrease the value in the vicinity of the coast 
    476                   presto(ji,jj,jk) = presto(ji,jj,1 ) * 0.5_wp * (  1._wp - COS( rpi*zdct(ji,jj,jk)/zinfl)  ) 
    477                   !   ... Vertical variation from zsdmp (sea surface) to zbdmp (bottom) 
    478                   presto(ji,jj,jk) = presto(ji,jj,jk) * (  zbdmp + (zsdmp-zbdmp) * EXP(-fsdept(ji,jj,jk)/pn_dep)  ) 
    479                END DO 
    480             END DO 
    481          END DO 
    482          ! 
    483       ENDIF 
    484  
    485       !                                  ! ========================= 
    486       !                                  !  Med and Red Sea damping    (ORCA configuration only) 
    487       !                                  ! ========================= 
    488       IF( cp_cfg == "orca" .AND. ( kn_hdmp > 0 .OR. kn_hdmp == -1 ) ) THEN 
    489          IF(lwp)WRITE(numout,*) 
    490          IF(lwp)WRITE(numout,*) '              ORCA configuration: Damping in Med and Red Seas' 
    491          ! 
    492          zmrs(:,:) = 0._wp 
    493          ! 
    494          SELECT CASE ( jp_cfg ) 
    495          !                                           ! ======================= 
    496          CASE ( 4 )                                  !  ORCA_R4 configuration  
    497             !                                        ! ======================= 
    498             ij0 =  50   ;   ij1 =  56                    ! Mediterranean Sea 
    499  
    500             ii0 =  81   ;   ii1 =  91   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    501             ij0 =  50   ;   ij1 =  55 
    502             ii0 =  75   ;   ii1 =  80   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    503             ij0 =  52   ;   ij1 =  53 
    504             ii0 =  70   ;   ii1 =  74   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    505             ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea 
    506             DO jk = 1, 17 
    507                zhfac (jk) = 0.5_wp * (  1._wp - COS( rpi * REAL(jk-1,wp) / 16._wp )  ) / rday 
    508             END DO 
    509             DO jk = 18, jpkm1 
    510                zhfac (jk) = 1._wp / rday 
    511             END DO 
    512             !                                        ! ======================= 
    513          CASE ( 2 )                                  !  ORCA_R2 configuration  
    514             !                                        ! ======================= 
    515             ij0 =  96   ;   ij1 = 110                    ! Mediterranean Sea 
    516             ii0 = 157   ;   ii1 = 181   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    517             ij0 = 100   ;   ij1 = 110 
    518             ii0 = 144   ;   ii1 = 156   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    519             ij0 = 100   ;   ij1 = 103 
    520             ii0 = 139   ;   ii1 = 143   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    521             ! 
    522             ij0 = 101   ;   ij1 = 102                    ! Decrease before Gibraltar Strait 
    523             ii0 = 139   ;   ii1 = 141   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0._wp 
    524             ii0 = 142   ;   ii1 = 142   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp 
    525             ii0 = 143   ;   ii1 = 143   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp 
    526             ii0 = 144   ;   ii1 = 144   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75_wp 
    527             ! 
    528             ij0 =  87   ;   ij1 =  96                    ! Red Sea 
    529             ii0 = 147   ;   ii1 = 163   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    530             ! 
    531             ij0 =  91   ;   ij1 =  91                    ! Decrease before Bab el Mandeb Strait 
    532             ii0 = 153   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.80_wp 
    533             ij0 =  90   ;   ij1 =  90 
    534             ii0 = 153   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp 
    535             ij0 =  89   ;   ij1 =  89 
    536             ii0 = 158   ;   ii1 = 160   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp 
    537             ij0 =  88   ;   ij1 =  88 
    538             ii0 = 160   ;   ii1 = 163   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0._wp 
    539             ! Smooth transition from 0 at surface to 1./rday at the 18th level in Med and Red Sea 
    540             DO jk = 1, 17 
    541                zhfac (jk) = 0.5_wp * (  1._wp - COS( rpi * REAL(jk-1,wp) / 16._wp )  ) / rday 
    542             END DO 
    543             DO jk = 18, jpkm1 
    544                zhfac (jk) = 1._wp / rday 
    545             END DO 
    546             !                                        ! ======================= 
    547          CASE ( 05 )                                 !  ORCA_R05 configuration 
    548             !                                        ! ======================= 
    549             ii0 = 568   ;   ii1 = 574                    ! Mediterranean Sea 
    550             ij0 = 324   ;   ij1 = 333   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    551             ii0 = 575   ;   ii1 = 658 
    552             ij0 = 314   ;   ij1 = 366   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    553             ! 
    554             ii0 = 641   ;   ii1 = 651                    ! Black Sea (remaining part 
    555             ij0 = 367   ;   ij1 = 372   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    556             ! 
    557             ij0 = 324   ;   ij1 = 333                    ! Decrease before Gibraltar Strait 
    558             ii0 = 565   ;   ii1 = 565   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp / 90._wp 
    559             ii0 = 566   ;   ii1 = 566   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40_wp 
    560             ii0 = 567   ;   ii1 = 567   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.75_wp 
    561             ! 
    562             ii0 = 641   ;   ii1 = 665                    ! Red Sea 
    563             ij0 = 270   ;   ij1 = 310   ;   zmrs( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1._wp 
    564             ! 
    565             ii0 = 666   ;   ii1 = 675                    ! Decrease before Bab el Mandeb Strait 
    566             ij0 = 270   ;   ij1 = 290    
    567             DO ji = mi0(ii0), mi1(ii1) 
    568                zmrs( ji , mj0(ij0):mj1(ij1) ) = 0.1_wp * ABS( FLOAT(ji - mi1(ii1)) ) 
    569             END DO  
    570             zsdmp = 1._wp / ( pn_surf * rday ) 
    571             zbdmp = 1._wp / ( pn_bot  * rday ) 
    572             DO jk = 1, jpk 
    573                zhfac(jk) = (  zbdmp + (zsdmp-zbdmp) * EXP( -fsdept(1,1,jk)/pn_dep )  ) 
    574             END DO 
    575             !                                       ! ======================== 
    576          CASE ( 025 )                               !  ORCA_R025 configuration  
    577             !                                       ! ======================== 
    578             CALL ctl_stop( ' Not yet implemented in ORCA_R025' ) 
    579             ! 
    580          END SELECT 
    581  
    582          DO jk = 1, jpkm1 
    583             presto(:,:,jk) = zmrs(:,:) * zhfac(jk) + ( 1._wp - zmrs(:,:) ) * presto(:,:,jk) 
    584          END DO 
    585  
    586          ! Mask resto array and set to 0 first and last levels 
    587          presto(:,:, : ) = presto(:,:,:) * tmask(:,:,:) 
    588          presto(:,:, 1 ) = 0._wp 
    589          presto(:,:,jpk) = 0._wp 
    590          !                         !--------------------! 
    591       ELSE                         !     No damping     ! 
    592          !                         !--------------------! 
    593          CALL ctl_stop( 'Choose a correct value of nn_hdmp or put ln_tradmp to FALSE' ) 
    594       ENDIF 
    595 #endif 
    596  
    597       !                            !--------------------------------! 
    598       IF( kn_file == 1 ) THEN      !  save damping coef. in a file  ! 
    599          !                         !--------------------------------! 
    600          IF(lwp) WRITE(numout,*) '              create damping.coeff.nc file' 
    601          IF( cdtype == 'TRA' ) cfile = 'damping.coeff' 
    602          IF( cdtype == 'TRC' ) cfile = 'damping.coeff.trc' 
    603          IF( cdtype == 'DYN' ) cfile = 'damping.coeff.dyn' 
    604          cfile = TRIM( cfile ) 
    605          CALL iom_open  ( cfile, inum0, ldwrt = .TRUE., kiolib = jprstlib ) 
    606          CALL iom_rstput( 0, 0, inum0, 'Resto', presto ) 
    607          CALL iom_close ( inum0 ) 
    608       ENDIF 
    609       ! 
    610       CALL wrk_dealloc( jpk, zhfac) 
    611       CALL wrk_dealloc( jpi, jpj, zmrs ) 
    612       CALL wrk_dealloc( jpi, jpj, jpk, zdct ) 
    613       ! 
    614       IF( nn_timing == 1 )  CALL timing_stop('dtacof') 
    615       ! 
    616    END SUBROUTINE dtacof 
    617  
    618  
    619    SUBROUTINE cofdis( pdct ) 
    620       !!---------------------------------------------------------------------- 
    621       !!                 ***  ROUTINE cofdis  *** 
    622       !! 
    623       !! ** Purpose :   Compute the distance between ocean T-points and the 
    624       !!      ocean model coastlines. Save the distance in a NetCDF file. 
    625       !! 
    626       !! ** Method  :   For each model level, the distance-to-coast is  
    627       !!      computed as follows :  
    628       !!       - The coastline is defined as the serie of U-,V-,F-points 
    629       !!      that are at the ocean-land bound. 
    630       !!       - For each ocean T-point, the distance-to-coast is then  
    631       !!      computed as the smallest distance (on the sphere) between the  
    632       !!      T-point and all the coastline points. 
    633       !!       - For land T-points, the distance-to-coast is set to zero. 
    634       !!      C A U T I O N : Computation not yet implemented in mpp case. 
    635       !! 
    636       !! ** Action  : - pdct, distance to the coastline (argument) 
    637       !!              - NetCDF file 'dist.coast.nc'  
    638       !!---------------------------------------------------------------------- 
    639       USE ioipsl      ! IOipsl librairy 
    640       !! 
    641       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) ::   pdct   ! distance to the coastline 
    642       !! 
    643       INTEGER ::   ji, jj, jk, jl   ! dummy loop indices 
    644       INTEGER ::   iju, ijt, icoast, itime, ierr, icot   ! local integers 
    645       CHARACTER (len=32) ::   clname                     ! local name 
    646       REAL(wp) ::   zdate0                               ! local scalar 
    647       REAL(wp), POINTER, DIMENSION(:,:) ::  zxt, zyt, zzt, zmask 
    648       REAL(wp), POINTER, DIMENSION(:  ) ::  zxc, zyc, zzc, zdis    ! temporary workspace 
    649       LOGICAL , ALLOCATABLE, DIMENSION(:,:) ::  llcotu, llcotv, llcotf   ! 2D logical workspace 
    650       !!---------------------------------------------------------------------- 
    651       ! 
    652       IF( nn_timing == 1 )  CALL timing_start('cofdis') 
    653       ! 
    654       CALL wrk_alloc( jpi, jpj , zxt, zyt, zzt, zmask    ) 
    655       CALL wrk_alloc( 3*jpi*jpj, zxc, zyc, zzc, zdis     ) 
    656       ALLOCATE( llcotu(jpi,jpj), llcotv(jpi,jpj), llcotf(jpi,jpj)  ) 
    657       ! 
    658       IF( lk_mpp    )   CALL mpp_sum( ierr ) 
    659       IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'cofdis: requested local arrays unavailable') 
    660  
    661       ! 0. Initialization 
    662       ! ----------------- 
    663       IF(lwp) WRITE(numout,*) 
    664       IF(lwp) WRITE(numout,*) 'cofdis : compute the distance to coastline' 
    665       IF(lwp) WRITE(numout,*) '~~~~~~' 
    666       IF(lwp) WRITE(numout,*) 
    667       IF( lk_mpp ) & 
    668            & CALL ctl_stop('         Computation not yet implemented with key_mpp_...', & 
    669            &               '         Rerun the code on another computer or ', & 
    670            &               '         create the "dist.coast.nc" file using IDL' ) 
    671  
    672       pdct(:,:,:) = 0._wp 
    673       zxt(:,:) = COS( rad * gphit(:,:) ) * COS( rad * glamt(:,:) ) 
    674       zyt(:,:) = COS( rad * gphit(:,:) ) * SIN( rad * glamt(:,:) ) 
    675       zzt(:,:) = SIN( rad * gphit(:,:) ) 
    676  
    677  
    678       ! 1. Loop on vertical levels 
    679       ! -------------------------- 
    680       !                                                ! =============== 
    681       DO jk = 1, jpkm1                                 ! Horizontal slab 
    682          !                                             ! =============== 
    683          ! Define the coastline points (U, V and F) 
    684          DO jj = 2, jpjm1 
    685             DO ji = 2, jpim1 
    686                zmask(ji,jj) =  ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 
    687                    &           + tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk) ) 
    688                llcotu(ji,jj) = ( tmask(ji,jj,  jk) + tmask(ji+1,jj  ,jk) == 1._wp )  
    689                llcotv(ji,jj) = ( tmask(ji,jj  ,jk) + tmask(ji  ,jj+1,jk) == 1._wp )  
    690                llcotf(ji,jj) = ( zmask(ji,jj) > 0._wp ) .AND. ( zmask(ji,jj) < 4._wp ) 
    691             END DO 
    692          END DO 
    693  
    694          ! Lateral boundaries conditions 
    695          llcotu(:, 1 ) = umask(:,  2  ,jk) == 1 
    696          llcotu(:,jpj) = umask(:,jpjm1,jk) == 1 
    697          llcotv(:, 1 ) = vmask(:,  2  ,jk) == 1 
    698          llcotv(:,jpj) = vmask(:,jpjm1,jk) == 1 
    699          llcotf(:, 1 ) = fmask(:,  2  ,jk) == 1 
    700          llcotf(:,jpj) = fmask(:,jpjm1,jk) == 1 
    701  
    702          IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 
    703             llcotu( 1 ,:) = llcotu(jpim1,:) 
    704             llcotu(jpi,:) = llcotu(  2  ,:) 
    705             llcotv( 1 ,:) = llcotv(jpim1,:) 
    706             llcotv(jpi,:) = llcotv(  2  ,:) 
    707             llcotf( 1 ,:) = llcotf(jpim1,:) 
    708             llcotf(jpi,:) = llcotf(  2  ,:) 
    709          ELSE 
    710             llcotu( 1 ,:) = umask(  2  ,:,jk) == 1 
    711             llcotu(jpi,:) = umask(jpim1,:,jk) == 1 
    712             llcotv( 1 ,:) = vmask(  2  ,:,jk) == 1 
    713             llcotv(jpi,:) = vmask(jpim1,:,jk) == 1 
    714             llcotf( 1 ,:) = fmask(  2  ,:,jk) == 1 
    715             llcotf(jpi,:) = fmask(jpim1,:,jk) == 1 
    716          ENDIF 
    717          IF( nperio == 3 .OR. nperio == 4 ) THEN 
    718             DO ji = 1, jpim1 
    719                iju = jpi - ji + 1 
    720                llcotu(ji,jpj  ) = llcotu(iju,jpj-2) 
    721                llcotf(ji,jpjm1) = llcotf(iju,jpj-2) 
    722                llcotf(ji,jpj  ) = llcotf(iju,jpj-3) 
    723             END DO 
    724             DO ji = jpi/2, jpim1 
    725                iju = jpi - ji + 1 
    726                llcotu(ji,jpjm1) = llcotu(iju,jpjm1) 
    727             END DO 
    728             DO ji = 2, jpi 
    729                ijt = jpi - ji + 2 
    730                llcotv(ji,jpjm1) = llcotv(ijt,jpj-2) 
    731                llcotv(ji,jpj  ) = llcotv(ijt,jpj-3) 
    732             END DO 
    733          ENDIF 
    734          IF( nperio == 5 .OR. nperio == 6 ) THEN 
    735             DO ji = 1, jpim1 
    736                iju = jpi - ji 
    737                llcotu(ji,jpj  ) = llcotu(iju,jpjm1) 
    738                llcotf(ji,jpj  ) = llcotf(iju,jpj-2) 
    739             END DO 
    740             DO ji = jpi/2, jpim1 
    741                iju = jpi - ji 
    742                llcotf(ji,jpjm1) = llcotf(iju,jpjm1) 
    743             END DO 
    744             DO ji = 1, jpi 
    745                ijt = jpi - ji + 1 
    746                llcotv(ji,jpj  ) = llcotv(ijt,jpjm1) 
    747             END DO 
    748             DO ji = jpi/2+1, jpi 
    749                ijt = jpi - ji + 1 
    750                llcotv(ji,jpjm1) = llcotv(ijt,jpjm1) 
    751             END DO 
    752          ENDIF 
    753  
    754          ! Compute cartesian coordinates of coastline points 
    755          ! and the number of coastline points 
    756          icoast = 0 
    757          DO jj = 1, jpj 
    758             DO ji = 1, jpi 
    759                IF( llcotf(ji,jj) ) THEN 
    760                   icoast = icoast + 1 
    761                   zxc(icoast) = COS( rad*gphif(ji,jj) ) * COS( rad*glamf(ji,jj) ) 
    762                   zyc(icoast) = COS( rad*gphif(ji,jj) ) * SIN( rad*glamf(ji,jj) ) 
    763                   zzc(icoast) = SIN( rad*gphif(ji,jj) ) 
    764                ENDIF 
    765                IF( llcotu(ji,jj) ) THEN 
    766                   icoast = icoast+1 
    767                   zxc(icoast) = COS( rad*gphiu(ji,jj) ) * COS( rad*glamu(ji,jj) ) 
    768                   zyc(icoast) = COS( rad*gphiu(ji,jj) ) * SIN( rad*glamu(ji,jj) ) 
    769                   zzc(icoast) = SIN( rad*gphiu(ji,jj) ) 
    770                ENDIF 
    771                IF( llcotv(ji,jj) ) THEN 
    772                   icoast = icoast+1 
    773                   zxc(icoast) = COS( rad*gphiv(ji,jj) ) * COS( rad*glamv(ji,jj) ) 
    774                   zyc(icoast) = COS( rad*gphiv(ji,jj) ) * SIN( rad*glamv(ji,jj) ) 
    775                   zzc(icoast) = SIN( rad*gphiv(ji,jj) ) 
    776                ENDIF 
    777             END DO 
    778          END DO 
    779  
    780          ! Distance for the T-points 
    781          DO jj = 1, jpj 
    782             DO ji = 1, jpi 
    783                IF( tmask(ji,jj,jk) == 0._wp ) THEN 
    784                   pdct(ji,jj,jk) = 0._wp 
    785                ELSE 
    786                   DO jl = 1, icoast 
    787                      zdis(jl) = ( zxt(ji,jj) - zxc(jl) )**2   & 
    788                         &     + ( zyt(ji,jj) - zyc(jl) )**2   & 
    789                         &     + ( zzt(ji,jj) - zzc(jl) )**2 
    790                   END DO 
    791                   pdct(ji,jj,jk) = ra * SQRT( MINVAL( zdis(1:icoast) ) ) 
    792                ENDIF 
    793             END DO 
    794          END DO 
    795          !                                                ! =============== 
    796       END DO                                              !   End of slab 
    797       !                                                   ! =============== 
    798  
    799  
    800       ! 2. Create the  distance to the coast file in NetCDF format 
    801       ! ----------------------------------------------------------     
    802       clname = 'dist.coast' 
    803       itime  = 0 
    804       CALL ymds2ju( 0     , 1       , 1     , 0._wp , zdate0 ) 
    805       CALL restini( 'NONE', jpi     , jpj   , glamt, gphit ,   & 
    806          &          jpk   , gdept_1d, clname, itime, zdate0,   & 
    807          &          rdt   , icot                         ) 
    808       CALL restput( icot, 'Tcoast', jpi, jpj, jpk, 0, pdct ) 
    809       CALL restclo( icot ) 
    810       ! 
    811       CALL wrk_dealloc( jpi, jpj , zxt, zyt, zzt, zmask    ) 
    812       CALL wrk_dealloc( 3*jpi*jpj, zxc, zyc, zzc, zdis     ) 
    813       DEALLOCATE( llcotu, llcotv, llcotf  ) 
    814       ! 
    815       IF( nn_timing == 1 )  CALL timing_stop('cofdis') 
    816       ! 
    817    END SUBROUTINE cofdis 
    818    !!====================================================================== 
    819249END MODULE tradmp 
Note: See TracChangeset for help on using the changeset viewer.