Changeset 4739


Ignore:
Timestamp:
2014-08-13T10:46:04+02:00 (6 years ago)
Author:
timgraham
Message:

Updated C1D/dyndmp.F90 and trcdmp.F90 to read restoration coefficient from a file.
Modified namelist_top_ref to match new options
Bug fixes to DMP_TOOLS tool and addition of custom.F90 to allow users to make modifications. Also changed to use working precision (wp) throughout.

Location:
branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM
Files:
1 added
11 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/CONFIG/SHARED/namelist_top_ref

    r4340 r4739  
    7676&namtrc_dmp    !   passive tracer newtonian damping    
    7777!----------------------------------------------------------------------- 
    78    nn_hdmp_tr  =   -1      !  horizontal shape =-1, damping in Med and Red Seas only 
    79                            !                   =XX, damping poleward of XX degrees (XX>0) 
    80                            !                      + F(distance-to-coast) + Red and Med Seas 
    8178   nn_zdmp_tr  =    1      !  vertical   shape =0    damping throughout the water column 
    8279                           !                   =1 no damping in the mixing layer (kz  criteria) 
    8380                           !                   =2 no damping in the mixed  layer (rho crieria) 
    84    rn_surf_tr  =   50.     !  surface time scale of damping   [days] 
    85    rn_bot_tr   =  360.     !  bottom  time scale of damping   [days] 
    86    rn_dep_tr   =  800.     !  depth of transition between rn_surf and rn_bot [meters] 
    87    nn_file_tr  =    0      !  create a damping.coeff NetCDF file (=1) or not (=0) 
     81   cn_resto  = 'resto_tr.nc'    !  create a damping.coeff NetCDF file (=1) or not (=0) 
    8882/ 
    8983!----------------------------------------------------------------------- 
  • branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/OPA_SRC/C1D/dyndmp.F90

    r4624 r4739  
    33   !!                       ***  MODULE  dyndmp  *** 
    44   !! Ocean dynamics: internal restoring trend on momentum (U and V current) 
     5   !!                 This should only be used for C1D case in current form 
    56   !!====================================================================== 
    67   !! History :  3.5  ! 2013-08  (D. Calvert)  Original code 
     8   !!            3.6  ! 2014-08  (T. Graham) Modified to use netcdf file of 
     9   !!                             restoration coefficients supplied to tradmp 
    710   !!---------------------------------------------------------------------- 
    811 
     
    2730   USE wrk_nemo       ! Memory allocation 
    2831   USE timing         ! Timing 
     32   USE iom            ! I/O manager 
    2933 
    3034   IMPLICIT NONE 
     
    7579      NAMELIST/namc1d_dyndmp/ ln_dyndmp 
    7680      INTEGER :: ios 
     81      INTEGER :: imask 
    7782      !!---------------------------------------------------------------------- 
    7883 
     
    9398         WRITE(numout,*) '      add a damping term or not       ln_dyndmp = ', ln_dyndmp 
    9499         WRITE(numout,*) '   Namelist namtra_dmp    : Set damping parameters' 
    95          WRITE(numout,*) '      horizontal damping option       nn_hdmp   = ', nn_hdmp 
    96          WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp, '(non-C1D zoom: forced to 0)' 
    97          WRITE(numout,*) '      surface time scale (days)       rn_surf   = ', rn_surf 
    98          WRITE(numout,*) '      bottom time scale (days)        rn_bot    = ', rn_bot 
    99          WRITE(numout,*) '      depth of transition (meters)    rn_dep    = ', rn_dep 
    100          WRITE(numout,*) '      create a damping.coeff file     nn_file   = ', nn_file 
     100         WRITE(numout,*) '      Apply relaxation   or not       ln_tradmp = ', ln_tradmp 
     101         WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp 
     102         WRITE(numout,*) '      Damping file name               cn_resto  = ', cn_resto 
    101103         WRITE(numout,*) 
    102104      ENDIF 
     
    106108         IF( dyn_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dyn_dmp_init: unable to allocate arrays' ) 
    107109         ! 
    108 #if ! defined key_c1d 
    109          SELECT CASE ( nn_hdmp )             !==   control print of horizontal option   ==! 
    110          CASE (  -1  )   ;   IF(lwp) WRITE(numout,*) '   momentum damping in the Med & Red seas only' 
    111          CASE ( 1:90 )   ;   IF(lwp) WRITE(numout,*) '   momentum damping poleward of', nn_hdmp, ' degrees' 
    112          CASE DEFAULT 
    113             WRITE(ctmp1,*) '          bad flag value for nn_hdmp = ', nn_hdmp 
    114             CALL ctl_stop(ctmp1) 
    115          END SELECT 
    116          ! 
    117 #endif 
    118110         SELECT CASE ( nn_zdmp )             !==   control print of vertical option   ==! 
    119111         CASE ( 0    )   ;   IF(lwp) WRITE(numout,*) '   momentum damping throughout the water column' 
     
    132124         utrdmp(:,:,:) = 0._wp               ! internal damping trends 
    133125         vtrdmp(:,:,:) = 0._wp 
    134          !                                   !==   Damping coefficients calculation:                           ==! 
    135          !                                   !==   use tradmp.F90 subroutines dtacof, dtacof_zoom and cofdis   ==! 
    136          !                                   !!!   NOTE: these need to be altered for use in this module if  
    137          !                                   !!!       they are to be used outside the C1D context  
    138          !                                   !!!       (use of U,V grid variables) 
    139          IF( lzoom .AND. .NOT. lk_c1d ) THEN   ;   CALL dtacof_zoom( resto_uv ) 
    140          ELSE               ;   CALL dtacof( nn_hdmp, rn_surf, rn_bot, rn_dep, nn_file, 'DYN', resto_uv ) 
    141          ENDIF 
    142          ! 
     126         ! 
     127         !Read in mask from file 
     128         CALL iom_open ( cn_resto, imask) 
     129         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto) 
     130         CALL iom_close( imask ) 
    143131      ENDIF 
    144132      ! 
  • branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90

    r4738 r4739  
    3636   USE wrk_nemo       ! Memory allocation 
    3737   USE timing         ! Timing 
     38   USE iom 
    3839 
    3940   IMPLICIT NONE 
     
    4647   LOGICAL , PUBLIC ::   ln_tradmp   !: internal damping flag 
    4748   INTEGER , PUBLIC ::   nn_zdmp     ! = 0/1/2 flag for damping in the mixed layer 
    48    CHARACTER(LEN=200) :: cn_resto      ! name of netcdf file containing restoration coefficient field 
     49   CHARACTER(LEN=200) , PUBLIC :: cn_resto      ! name of netcdf file containing restoration coefficient field 
    4950   LOGICAL , PUBLIC ::   ln_miss      ! check for missing data in T/S data file (slow?) 
    5051   REAL(wp), PUBLIC ::   rn_miss      ! Value of missing data 
     52   ! 
    5153 
    5254 
     
    189191      !!---------------------------------------------------------------------- 
    190192 
    191       NAMELIST/namtra_dmp/ ln_tradmp, rn_dmp, nn_zdmp, cn_resto 
     193      NAMELIST/namtra_dmp/ ln_tradmp, nn_zdmp, cn_resto 
    192194      INTEGER ::  ios         ! Local integer for output status of namelist read 
    193       REAL(wp), POINTER, DIMENSION(:,:,:) ::  dmp_mask   ! 3D mask 
     195      INTEGER :: imask        ! File handle  
    194196      !!---------------------------------------------------------------------- 
    195197 
     
    224226         CASE ( 1 )  ; IF(lwp) WRITE(numout,*) '   no tracer damping in the turbocline' 
    225227         CASE ( 2 )  ; IF(lwp) WRITE(numout,*) '   no tracer damping in the mixed layer' 
     228         END SELECT 
    226229 
    227230         !Initialisation of dtatsd - Would it be better to have dmpdta routine 
  • branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90

    r4359 r4739  
    302302      !!---------------------------------------------------------------------- 
    303303      ! 
     304      INTEGER, imask  !local file handle 
     305 
    304306      IF( nn_timing == 1 )  CALL timing_start('trc_dmp_init') 
    305307      ! 
    306       SELECT CASE ( nn_hdmp_tr ) 
    307       CASE (  -1  )   ;   IF(lwp) WRITE(numout,*) '   tracer damping in the Med & Red seas only' 
    308       CASE ( 1:90 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping poleward of', nn_hdmp_tr, ' degrees' 
    309       CASE DEFAULT 
    310          WRITE(ctmp1,*) '          bad flag value for nn_hdmp_tr = ', nn_hdmp_tr 
    311          CALL ctl_stop(ctmp1) 
    312       END SELECT 
    313308 
    314309      IF( lzoom )   nn_zdmp_tr = 0           ! restoring to climatology at closed north or south boundaries 
     
    325320         &   CALL ctl_stop( 'passive trace damping need key_tradmp to compute damping coef.' ) 
    326321      ! 
    327       !                          ! Damping coefficients initialization 
    328       IF( lzoom ) THEN   ;   CALL dtacof_zoom( restotr ) 
    329       ELSE               ;   CALL dtacof( nn_hdmp_tr, rn_surf_tr, rn_bot_tr, rn_dep_tr,  & 
    330                              &            nn_file_tr, 'TRC'     , restotr                ) 
    331       ENDIF 
     322      !                          ! Read damping coefficients from file 
     323      !Read in mask from file 
     324      CALL iom_open ( cn_resto_tr, imask) 
     325      CALL iom_get  ( imask, jpdom_autoglo, 'resto', restotr) 
     326      CALL iom_close( imask ) 
    332327      ! 
    333328      IF( nn_timing == 1 )  CALL timing_stop('trc_dmp_init') 
  • branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/NEMO/TOP_SRC/TRP/trcnam_trp.F90

    r4624 r4739  
    5151   !                                                 !!: ** newtonian damping namelist (nam_trcdmp) ** 
    5252   !                          !!* Namelist namtrc_dmp : passive tracer newtonian damping * 
    53    INTEGER , PUBLIC ::   nn_hdmp_tr    ! = 0/-1/'latitude' for damping over passive tracer 
    5453   INTEGER , PUBLIC ::   nn_zdmp_tr    ! = 0/1/2 flag for damping in the mixed layer 
    55    REAL(wp), PUBLIC ::   rn_surf_tr    ! surface time scale for internal damping        [days] 
    56    REAL(wp), PUBLIC ::   rn_bot_tr     ! bottom time scale for internal damping         [days] 
    57    REAL(wp), PUBLIC ::   rn_dep_tr     ! depth of transition between rn_surf and rn_bot [meters] 
    58    INTEGER , PUBLIC ::   nn_file_tr    ! = 1 create a damping.coeff NetCDF file 
     54   CHARACTER(LEN=200) , PUBLIC :: cn_resto_tr    !File containing restoration coefficient 
    5955 
    6056   !!---------------------------------------------------------------------- 
     
    8278      NAMELIST/namtrc_zdf/ ln_trczdf_exp  , nn_trczdf_exp 
    8379      NAMELIST/namtrc_rad/ ln_trcrad 
    84       NAMELIST/namtrc_dmp/ nn_hdmp_tr, nn_zdmp_tr, rn_surf_tr, & 
    85         &                  rn_bot_tr , rn_dep_tr , nn_file_tr 
     80      NAMELIST/namtrc_dmp/ nn_zdmp_tr , nn_file_tr 
    8681      !!---------------------------------------------------------------------- 
    8782 
     
    184179         WRITE(numout,*) '~~~~~~~' 
    185180         WRITE(numout,*) '   Namelist namtrc_dmp : set damping parameter' 
    186          WRITE(numout,*) '      tracer damping option          nn_hdmp_tr = ', nn_hdmp_tr 
    187181         WRITE(numout,*) '      mixed layer damping option     nn_zdmp_tr = ', nn_zdmp_tr, '(zoom: forced to 0)' 
    188          WRITE(numout,*) '      surface time scale (days)      rn_surf_tr = ', rn_surf_tr 
    189          WRITE(numout,*) '      bottom time scale (days)       rn_bot_tr  = ', rn_bot_tr 
    190          WRITE(numout,*) '      depth of transition (meters)   rn_dep_tr  = ', rn_dep_tr 
    191          WRITE(numout,*) '      create a damping.coeff file    nn_file_tr = ', nn_file_tr 
     182         WRITE(numout,*) '      Restoration coeff file    cn_resto_tr = ', cn_resto_tr 
    192183      ENDIF 
    193184      ! 
  • branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/namelist

    r4738 r4739  
    33    cp_cfz = 'antarctic'       ! Name of zoom configuration (arctic and antarctic have some special treatment if lzoom=.true.) 
    44    jp_cfg = 2                 ! Resolution of the model (used for med_red_seas damping) 
    5     lzoom = .true.             ! Zoom configuration or not 
     5    lzoom = .false.            ! Zoom configuration or not 
    66    ln_full_field = .false.    ! Calculate coefficient over whole of domain  
    77    ln_med_red_seas = .true.   ! Damping in Med/Red Seas (or local modifications here if ln_full_field=.true.) 
     8    ln_old_31_lev_code = .true.   ! Replicate behaviour of old online code for 31 level model (Med/Red seas damping based on level number instead of depth) 
    89    ln_coast = .true.          ! Reduce near to coastlines  
    910    ln_zero_top_layer = .true. ! No damping in top layer 
  • branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/coast_dist.F90

    r4738 r4739  
    1919      !! 
    2020      IMPLICIT NONE 
    21       REAL(8), DIMENSION(jpi,jpj), INTENT( inout ) :: presto 
    22       REAL(8), DIMENSION(jpi,jpj) :: zdct 
    23       REAL(8) :: zinfl = 1000.e3  ! Distance of influence of coast line (could be 
     21      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) :: presto 
     22      REAL(wp), DIMENSION(jpi,jpj) :: zdct 
     23      REAL(wp) :: zinfl = 1000.e3_wp  ! Distance of influence of coast line (could be 
    2424                                  ! a namelist setting) 
    2525      INTEGER :: jj, ji           ! dummy loop indices 
     
    3030         DO ji = 1, jpi 
    3131            zdct(ji,jj) = MIN( zinfl, zdct(ji,jj) ) 
    32             presto(ji,jj) = presto(ji, jj) * 0.5 * (  1. - COS( rpi*zdct(ji,jj)/zinfl) ) 
     32            presto(ji,jj) = presto(ji, jj) * 0.5_wp * (  1._wp - COS( rpi*zdct(ji,jj)/zinfl) ) 
    3333         END DO 
    3434      END DO 
     
    5757      !!---------------------------------------------------------------------- 
    5858      !! 
    59       REAL(8), DIMENSION(jpi,jpj), INTENT( out ) ::   pdct   ! distance to the coastline 
     59      REAL(wp), DIMENSION(jpi,jpj), INTENT( out ) ::   pdct   ! distance to the coastline 
    6060      !! 
    6161      INTEGER ::   ji, jj, jl   ! dummy loop indices 
    6262      INTEGER ::   iju, ijt, icoast, itime, ierr, icot   ! local integers 
    6363      CHARACTER (len=32) ::   clname                     ! local name 
    64       REAL(8) ::   zdate0                               ! local scalar 
    65       REAL(8), POINTER, DIMENSION(:,:) ::  zxt, zyt, zzt, zmask 
    66       REAL(8), POINTER, DIMENSION(:  ) ::  zxc, zyc, zzc, zdis    ! temporary workspace 
     64      REAL(wp) ::   zdate0                               ! local scalar 
     65      REAL(wp), POINTER, DIMENSION(:,:) ::  zxt, zyt, zzt, zmask 
     66      REAL(wp), POINTER, DIMENSION(:  ) ::  zxc, zyc, zzc, zdis    ! temporary workspace 
    6767      LOGICAL , ALLOCATABLE, DIMENSION(:,:) ::  llcotu, llcotv, llcotf   ! 2D logical workspace 
    6868 
     
    9090      CALL check_nf90( nf90_get_var( ncin, fmask_id, fmask, (/ 1,1 /), (/ jpi, jpj /) ) ) 
    9191 
    92       pdct(:,:) = 0. 
     92      pdct(:,:) = 0._wp 
    9393      zxt(:,:) = COS( rad * gphit(:,:) ) * COS( rad * glamt(:,:) ) 
    9494      zyt(:,:) = COS( rad * gphit(:,:) ) * SIN( rad * glamt(:,:) ) 
     
    101101               zmask(ji,jj) =  ( tmask(ji,jj+1) + tmask(ji+1,jj+1) & 
    102102                   &           + tmask(ji,jj  ) + tmask(ji+1,jj  ) ) 
    103                llcotu(ji,jj) = ( tmask(ji,jj ) + tmask(ji+1,jj  ) == 1. )  
    104                llcotv(ji,jj) = ( tmask(ji,jj  ) + tmask(ji  ,jj+1) == 1. )  
    105                llcotf(ji,jj) = ( zmask(ji,jj) > 0. ) .AND. ( zmask(ji,jj) < 4. ) 
     103               llcotu(ji,jj) = ( tmask(ji,jj ) + tmask(ji+1,jj  ) == 1._wp )  
     104               llcotv(ji,jj) = ( tmask(ji,jj  ) + tmask(ji  ,jj+1) == 1._wp )  
     105               llcotf(ji,jj) = ( zmask(ji,jj) > 0._wp ) .AND. ( zmask(ji,jj) < 4._wp ) 
    106106            END DO 
    107107         END DO 
     
    196196         DO jj = 1, jpj 
    197197            DO ji = 1, jpi 
    198                IF( tmask(ji,jj) == 0. ) THEN 
    199                   pdct(ji,jj) = 0. 
     198               IF( tmask(ji,jj) == 0._wp ) THEN 
     199                  pdct(ji,jj) = 0._wp 
    200200               ELSE 
    201201                  DO jl = 1, icoast 
  • branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/make_dmp_file.F90

    r4738 r4739  
    2424  USE med_red_seas 
    2525  USE zoom 
     26  USE custom 
    2627 
    2728  IMPLICIT NONE 
    2829  INTEGER  :: ji, jj, jk                         ! dummpy loop variables 
    29   REAL(8) :: zsdmp, zbdmp                     ! Surface and bottom damping coeff 
     30  REAL(wp) :: zsdmp, zbdmp                     ! Surface and bottom damping coeff 
    3031  CHARACTER(LEN=200) :: meshfile = 'mesh_mask.nc'   ! mesh file 
    31   CHARACTER(LEN=200) :: outfile = 'dmp_mask.nc'     ! output file 
    32   REAL(8) :: zlat, zlat2, zlat0 
     32  CHARACTER(LEN=200) :: outfile = 'resto.nc'     ! output file 
     33  REAL(wp) :: zlat, zlat2, zlat0 
    3334 
    3435  ! Read namelist 
     
    5556 
    5657  !Calculate surface and bottom damping coefficients 
    57   zsdmp = 1. / ( pn_surf * rday ) 
    58   zbdmp = 1. / ( pn_bot  * rday ) 
     58  zsdmp = 1._wp / ( pn_surf * rday ) 
     59  zbdmp = 1._wp / ( pn_bot  * rday ) 
    5960 
    6061  !Loop through levels and read in tmask for each level as starting point for 
    6162  !coefficient array 
    6263  DO jk = 1, jpk-1 
    63      resto(:,:) = 0. 
     64     resto(:,:) = 0._wp 
    6465      
    6566     IF (.NOT. (jk == 1 .AND. ln_zero_top_layer) ) THEN  
     
    8384                    zlat = ABS(gphit(ji,jj)) 
    8485                    IF ( nn_hdmp <= zlat .AND. zlat <= zlat2 ) THEN 
    85                        resto(ji,jj) = resto(ji,jj) * 0.5 * (  1. - COS( rpi*(zlat-nn_hdmp)/zlat0 ) ) 
     86                       resto(ji,jj) = resto(ji,jj) * 0.5_wp * (  1._wp - COS( rpi*(zlat-nn_hdmp)/zlat0 ) ) 
    8687                    ELSE IF ( zlat < nn_hdmp ) THEN 
    87                        resto(ji,jj) = 0. 
     88                       resto(ji,jj) = 0._wp 
    8889                    ENDIF 
    8990                 END DO 
     
    9899 
    99100        ! Damping in Med/Red Seas (or local modifications if full field is set) 
    100         IF (ln_med_red_seas .AND. (cp_cfg == 'orca')) THEN 
    101            CALL med_red_dmp(resto) 
     101        IF (ln_med_red_seas .AND. (cp_cfg == 'orca') .AND. (.NOT. lzoom)) THEN 
     102           CALL med_red_dmp(resto, jk, ln_old_31_lev_code) 
    102103        ENDIF 
    103104 
     
    105106              CALL dtacof_zoom(resto, tmask) 
    106107        ENDIF 
    107  
     108         
     109        !Any user modifications can be added in the custom module 
     110        IF ( ln_custom ) THEN 
     111              CALL custom_resto( resto ) 
     112        ENDIF 
    108113     ENDIF 
    109114 
  • branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/med_red_seas.F90

    r4738 r4739  
    88   CONTAINS   
    99 
    10    SUBROUTINE med_red_dmp(presto) 
     10   SUBROUTINE med_red_dmp(presto, jk, ln_31_lev) 
    1111      !!------------------------------------ 
    1212      !!    **ROUTINE: med_red_dmp 
     
    1717      !!----------------------------------- 
    1818      INTEGER :: ij0,ij1,ii0,ii1,ji,jj      
    19       REAL(8), DIMENSION(:,:), ALLOCATABLE :: zmrs 
    20       REAL(8) :: zhfac, zsdmp, zbdmp 
    21       REAL(8), DIMENSION(jpi,jpj), INTENT(inout) :: presto 
     19      INTEGER, INTENT(in) :: jk 
     20      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zmrs 
     21      REAL(wp) :: zhfac, zsdmp, zbdmp 
     22      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: presto 
     23      LOGICAL, INTENT(in), OPTIONAL :: ln_31_lev 
     24      LOGICAL :: l_31_lev 
    2225 
    2326      WRITE(numout,*) 'ORCA Med and Red Seas Damping' 
    2427       
     28      IF ( PRESENT(ln_31_lev)) THEN 
     29         l_31_lev = ln_31_lev 
     30      ELSE 
     31         l_31_lev = .false. 
     32      ENDIF 
     33       
    2534      ALLOCATE( zmrs(jpi, jpj) ) 
    2635         ! 
    27          zmrs(:,:) = 0. 
     36         zmrs(:,:) = 0._wp 
    2837         ! 
    2938         SELECT CASE ( jp_cfg ) 
     
    4352            !                                        ! ======================= 
    4453            ij0 =  96   ;   ij1 = 110                    ! Mediterranean Sea 
    45             ii0 = 157   ;   ii1 = 181   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1. 
     54            ii0 = 157   ;   ii1 = 181   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp 
    4655            ij0 = 100   ;   ij1 = 110 
    47             ii0 = 144   ;   ii1 = 156   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1. 
     56            ii0 = 144   ;   ii1 = 156   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp 
    4857            ij0 = 100   ;   ij1 = 103 
    49             ii0 = 139   ;   ii1 = 143   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1. 
     58            ii0 = 139   ;   ii1 = 143   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp 
    5059            ! 
    5160            ij0 = 101   ;   ij1 = 102                    ! Decrease before Gibraltar Strait 
    52             ii0 = 139   ;   ii1 = 141   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0. 
    53             ii0 = 142   ;   ii1 = 142   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1. / 90. 
    54             ii0 = 143   ;   ii1 = 143   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.40 
    55             ii0 = 144   ;   ii1 = 144   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.75 
     61            ii0 = 139   ;   ii1 = 141   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0._wp 
     62            ii0 = 142   ;   ii1 = 142   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp / 90._wp 
     63            ii0 = 143   ;   ii1 = 143   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.40_wp 
     64            ii0 = 144   ;   ii1 = 144   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.75_wp 
    5665            ! 
    5766            ij0 =  87   ;   ij1 =  96                    ! Red Sea 
    58             ii0 = 147   ;   ii1 = 163   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1. 
     67            ii0 = 147   ;   ii1 = 163   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp 
    5968            ! 
    6069            ij0 =  91   ;   ij1 =  91                    ! Decrease before Bab el Mandeb Strait 
    61             ii0 = 153   ;   ii1 = 160   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.80 
     70            ii0 = 153   ;   ii1 = 160   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.80_wp 
    6271            ij0 =  90   ;   ij1 =  90 
    63             ii0 = 153   ;   ii1 = 160   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.40 
     72            ii0 = 153   ;   ii1 = 160   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.40_wp 
    6473            ij0 =  89   ;   ij1 =  89 
    65             ii0 = 158   ;   ii1 = 160   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1. / 90. 
     74            ii0 = 158   ;   ii1 = 160   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp / 90._wp 
    6675            ij0 =  88   ;   ij1 =  88 
    67             ii0 = 160   ;   ii1 = 163   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0. 
     76            ii0 = 160   ;   ii1 = 163   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0._wp 
    6877            ! 
    6978            !                                        ! ======================= 
     
    7180            !                                        ! ======================= 
    7281            ii0 = 568   ;   ii1 = 574                    ! Mediterranean Sea 
    73             ij0 = 324   ;   ij1 = 333   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1. 
     82            ij0 = 324   ;   ij1 = 333   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp 
    7483            ii0 = 575   ;   ii1 = 658 
    75             ij0 = 314   ;   ij1 = 366   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1. 
     84            ij0 = 314   ;   ij1 = 366   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp 
    7685            ! 
    7786            ii0 = 641   ;   ii1 = 651                    ! Black Sea (remaining part 
    78             ij0 = 367   ;   ij1 = 372   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1. 
     87            ij0 = 367   ;   ij1 = 372   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp 
    7988            ! 
    8089            ij0 = 324   ;   ij1 = 333                    ! Decrease before Gibraltar Strait 
    81             ii0 = 565   ;   ii1 = 565   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1. / 90. 
    82             ii0 = 566   ;   ii1 = 566   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.40 
    83             ii0 = 567   ;   ii1 = 567   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.75 
     90            ii0 = 565   ;   ii1 = 565   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp / 90._wp 
     91            ii0 = 566   ;   ii1 = 566   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.40_wp 
     92            ii0 = 567   ;   ii1 = 567   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 0.75_wp 
    8493            ! 
    8594            ii0 = 641   ;   ii1 = 665                    ! Red Sea 
    86             ij0 = 270   ;   ij1 = 310   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1. 
     95            ij0 = 270   ;   ij1 = 310   ;   zmrs( ii0:ii1 , ij0:ij1 ) = 1._wp 
    8796            ! 
    8897            ii0 = 666   ;   ii1 = 675                    ! Decrease before Bab el Mandeb Strait 
    8998            ij0 = 270   ;   ij1 = 290    
    9099            DO ji = ii0, ii1 
    91                zmrs( ji , ij0:ij1 ) = 0.1 * ABS( FLOAT(ji - ii1) ) 
     100               zmrs( ji , ij0:ij1 ) = 0.1_wp * ABS( FLOAT(ji - ii1) ) 
    92101            END DO  
    93102            !                                       ! ======================== 
     
    100109         END SELECT 
    101110 
    102          ! Note that the original "online" code had a dependency on model levels 
    103          ! here (as opposed to depth) 
    104          ! This has been removed but can be reproduced using the "custom" module 
    105          ! if required 
    106          zsdmp = 1. / ( pn_surf * rday ) 
    107          zbdmp = 1. / ( pn_bot  * rday ) 
    108          zhfac = (  zbdmp + (zsdmp-zbdmp) * EXP( -gdept(1,1)/pn_dep )  ) 
     111         zsdmp = 1._wp / ( pn_surf * rday ) 
     112         zbdmp = 1._wp / ( pn_bot  * rday ) 
    109113 
    110          presto(:,:) = zmrs(:,:) * zhfac + ( 1. - zmrs(:,:) ) * presto(:,:) 
     114         ! The l_31_lev option is used to reproduce the old behaviour of 
     115         ! defining the restoration coefficient based on the level number. 
     116         ! This is included to allow damping coefficients for reference 
     117         ! configurations to be kept the same. 
     118         IF (l_31_lev) THEN 
     119            IF (jk <= 17) THEN 
     120               zhfac = 0.5_wp * (  1. - COS( rpi * REAL(jk-1,wp) / 16._wp )  ) / rday 
     121            ELSE 
     122               zhfac = 1._wp / rday 
     123            ENDIF 
     124         ELSE 
     125            zhfac = (  zbdmp + (zsdmp-zbdmp) * EXP( -gdept(1,1)/pn_dep )  ) 
     126         ENDIF 
     127 
     128         presto(:,:) = zmrs(:,:) * zhfac + ( 1._wp - zmrs(:,:) ) * presto(:,:) 
    111129 
    112130         DEALLOCATE( zmrs )          
  • branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/utils.F90

    r4738 r4739  
    66  PUBLIC 
    77   
     8  INTEGER, PUBLIC, PARAMETER   :: dp=8 , sp=4, wp=dp 
    89  INTEGER :: tmask_id, umask_id, vmask_id, fmask_id 
    910  INTEGER :: gdept_id 
     
    1314  INTEGER :: jpi, jpj, jpk                      ! Size of domain 
    1415  INTEGER  :: ncin, ncout                              ! File handles for netCDF files 
    15   REAL(8), DIMENSION(:,:),   ALLOCATABLE :: gphit, glamt 
    16   REAL(8), DIMENSION(:,:),   ALLOCATABLE :: gphiu, glamu 
    17   REAL(8), DIMENSION(:,:),   ALLOCATABLE :: gphiv, glamv 
    18   REAL(8), DIMENSION(:,:),   ALLOCATABLE :: gphif, glamf 
    19   REAL(8), DIMENSION(:,:),   ALLOCATABLE :: tmask, umask, vmask, fmask 
    20   REAL(8), DIMENSION(:,:),   ALLOCATABLE :: gdept 
    21   REAL(8), DIMENSION(:,:),   ALLOCATABLE :: resto 
     16  REAL(wp), DIMENSION(:,:),   ALLOCATABLE :: gphit, glamt 
     17  REAL(wp), DIMENSION(:,:),   ALLOCATABLE :: gphiu, glamu 
     18  REAL(wp), DIMENSION(:,:),   ALLOCATABLE :: gphiv, glamv 
     19  REAL(wp), DIMENSION(:,:),   ALLOCATABLE :: gphif, glamf 
     20  REAL(wp), DIMENSION(:,:),   ALLOCATABLE :: tmask, umask, vmask, fmask 
     21  REAL(wp), DIMENSION(:,:),   ALLOCATABLE :: gdept 
     22  REAL(wp), DIMENSION(:,:),   ALLOCATABLE :: resto 
    2223 
    2324  INTEGER,PARAMETER :: numout = 6 
    2425  INTEGER,PARAMETER :: numerr = 0 
    2526  INTEGER,PARAMETER :: numnam = 11 
    26   REAL(8),PARAMETER :: rday = 86400           ! seconds in a day 
    27   REAL(8),PARAMETER ::   rpi = 3.141592653589793 
    28   REAL(8),PARAMETER ::   rad = 3.141592653589793/180. 
    29   REAL(8),PARAMETER ::   ra =  6371229. 
     27  REAL(wp),PARAMETER :: rday = 86400           ! seconds in a day 
     28  REAL(wp),PARAMETER ::   rpi = 3.141592653589793 
     29  REAL(wp),PARAMETER ::   rad = 3.141592653589793/180. 
     30  REAL(wp),PARAMETER ::   ra =  6371229. 
    3031 
    3132  ! Namelist variables 
     
    4243  LOGICAL :: ln_full_field = .true. 
    4344  LOGICAL :: ln_med_red_seas = .false. 
     45  LOGICAL :: ln_old_31_lev_code = .false. 
    4446  LOGICAL :: ln_zero_top_layer = .false. 
    4547  LOGICAL :: ln_custom = .false. 
    4648 
    4749  NAMELIST/nam_dmp_create/cp_cfg, cp_cfz, jp_cfg, lzoom, ln_full_field, & 
    48                           ln_med_red_seas, ln_coast, ln_zero_top_layer, ln_custom, & 
     50                          ln_med_red_seas, ln_old_31_lev_code, ln_coast, & 
     51                          ln_zero_top_layer, ln_custom, & 
    4952                          pn_surf, pn_bot, pn_dep, nn_hdmp, jperio 
    5053 
     
    106109     CALL check_nf90( nf90_def_dim(ncout, 'z', jpk, id_z) ) 
    107110 
    108      CALL check_nf90( nf90_def_var(ncout, 'resto', nf90_float, (/id_x,id_y,id_z/), resto_id ) ) 
     111     CALL check_nf90( nf90_def_var(ncout, 'resto', nf90_double, (/id_x,id_y,id_z/), resto_id ) ) 
    109112     CALL check_nf90( nf90_enddef(ncout) ) 
    110113 
  • branches/2014/dev_r4650_UKMO3_masked_damping/NEMOGCM/TOOLS/DMP_TOOLS/src/zoom.F90

    r4738 r4739  
    1818      !! ** Action  : - resto, the damping coeff. for T and S 
    1919      !!---------------------------------------------------------------------- 
    20       REAL(8), DIMENSION(jpi,jpj), INTENT(inout)  ::   presto   ! restoring coeff. (s-1) 
    21       REAL(8), DIMENSION(jpi,jpj), INTENT(in)  ::   mask   ! restoring coeff. (s-1) 
     20      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)  ::   presto   ! restoring coeff. (s-1) 
     21      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  ::   mask   ! restoring coeff. (s-1) 
    2222      ! 
    2323      INTEGER  ::   ji, jj, jn   ! dummy loop indices 
    24       REAL(8) ::   zlat, zlat0, zlat1, zlat2, z1_5d   ! local scalar 
    25       REAL(8), DIMENSION(6)  ::   zfact               ! 1Dworkspace 
     24      REAL(wp) ::   zlat, zlat0, zlat1, zlat2, z1_5d   ! local scalar 
     25      REAL(wp), DIMENSION(6)  ::   zfact               ! 1Dworkspace 
     26 
     27      !Namelist variables 
     28      LOGICAL :: lzoom_w, lzoom_e, lzoom_n, lzoom_s  
     29      NAMELIST/nam_zoom_dmp/lzoom_n,lzoom_e,lzoom_w,lzoom_s 
    2630      !!---------------------------------------------------------------------- 
    2731      ! 
     
    2933      ! 
    3034       
     35      ! Read namelist 
     36      OPEN( UNIT=numnam, FILE='namelist', FORM='FORMATTED', STATUS='OLD' ) 
     37      READ( numnam, nam_dmp_create ) 
     38      CLOSE( numnam ) 
    3139 
    32       zfact(1) =  1. 
    33       zfact(2) =  1. 
    34       zfact(3) = 11. / 12. 
    35       zfact(4) =  8. / 12. 
    36       zfact(5) =  4. / 12. 
    37       zfact(6) =  1. / 12. 
    38       zfact(:) = zfact(:) / ( 5. * rday )    ! 5 days max restoring time scale 
     40      zfact(1) =  1._wp 
     41      zfact(2) =  1._wp 
     42      zfact(3) = 11._wp / 12._wp 
     43      zfact(4) =  8._wp / 12._wp 
     44      zfact(5) =  4._wp / 12._wp 
     45      zfact(6) =  1._wp / 12._wp 
     46      zfact(:) = zfact(:) / ( 5._wp * rday )    ! 5 days max restoring time scale 
    3947 
    40       presto(:,:) = 0. 
     48      presto(:,:) = 0._wp 
    4149 
    4250      ! damping along the forced closed boundary over 6 grid-points 
     
    5159      IF( cp_cfz == "arctic" .OR. cp_cfz == "antarctic" ) THEN   !  ORCA configuration : arctic or antarctic zoom 
    5260         !                                        ! ==================================================== 
    53          IF(lwp) WRITE(numout,*) 
    54          IF(lwp .AND. cp_cfz == "arctic" ) WRITE(numout,*) '              dtacof_zoom : ORCA    Arctic zoom' 
    55          IF(lwp .AND. cp_cfz == "antarctic" ) WRITE(numout,*) '           dtacof_zoom : ORCA Antarctic zoom' 
    56          IF(lwp) WRITE(numout,*) 
     61         WRITE(numout,*) 
     62         IF(cp_cfz == "arctic" ) WRITE(numout,*) '              dtacof_zoom : ORCA    Arctic zoom' 
     63         IF(cp_cfz == "antarctic" ) WRITE(numout,*) '           dtacof_zoom : ORCA Antarctic zoom' 
     64         WRITE(numout,*) 
    5765         ! 
    5866         !                          ! Initialization :  
    59          presto(:,:) = 0. 
    60          zlat0 = 10.                     ! zlat0 : latitude strip where resto decreases 
    61          zlat1 = 30.                     ! zlat1 : resto = 1 before zlat1 
     67         presto(:,:) = 0._wp 
     68         zlat0 = 10._wp                     ! zlat0 : latitude strip where resto decreases 
     69         zlat1 = 30._wp                     ! zlat1 : resto = 1 before zlat1 
    6270         zlat2 = zlat1 + zlat0              ! zlat2 : resto decreases from 1 to 0 between zlat1 and zlat2 
    63          z1_5d = 1. / ( 5. * rday )   ! z1_5d : 1 / 5days 
     71         z1_5d = 1._wp / ( 5._wp * rday )   ! z1_5d : 1 / 5days 
    6472 
    6573         DO jj = 1, jpj 
     
    6775               zlat = ABS( gphit(ji,jj) ) 
    6876               IF( zlat1 <= zlat .AND. zlat <= zlat2 ) THEN 
    69                   presto(ji,jj) = 0.5 * z1_5d * (  1. - COS( rpi*(zlat2-zlat)/zlat0 )  )  
     77                  presto(ji,jj) = 0.5_wp * z1_5d * (  1._wp - COS( rpi*(zlat2-zlat)/zlat0 )  )  
    7078               ELSEIF( zlat < zlat1 ) THEN 
    7179                  presto(ji,jj) = z1_5d 
Note: See TracChangeset for help on using the changeset viewer.