Changeset 4738


Ignore:
Timestamp:
2014-08-11T17:47:48+02:00 (6 years ago)
Author:
timgraham
Message:

Modified tra_dmp module to read in restoration coefficient from a netcdf file

Added a tool to create the netcdf file - this replaces all of the hard coded resolution dependencies in tra_dmp_init

Location:
branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM
Files:
8 added
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/ARCH/OLD/arch-ifort_linux.fcm

    r3922 r4738  
    1717%NCDF_INC            -I/usr/local/netcdf/include  
    1818%NCDF_LIB            -L /usr/local/netcdf/lib -lnetcdf    
    19 %FC                  ifort 
     19%FC                  ifort  
    2020%FCFLAGS         -r8 -O3  -traceback  
    2121%FFLAGS       -r8 -O3  -traceback  
  • branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90

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