New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 5240 for branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC – NEMO

Ignore:
Timestamp:
2015-04-29T12:17:12+02:00 (9 years ago)
Author:
davestorkey
Message:

Update UKMO nn_etau_revision branch with trunk changes to rev 5107.

Location:
branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90

    r5239 r5240  
    125125            IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo 
    126126            IF(lwp) THEN  
    127                     WRITE(numout,*) '             Tidal cpt name    -     Phase speed (deg/hr)'             
     127                    WRITE(numout,*) '             Tidal components: '  
    128128               DO itide = 1, nb_harmo 
    129                   WRITE(numout,*)  '             ', Wave(ntide(itide))%cname_tide, omega_tide(itide)   
     129                  WRITE(numout,*)  '                 ', Wave(ntide(itide))%cname_tide  
    130130               END DO 
    131131            ENDIF  

  • branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/C1D/dyndmp.F90

    r5239 r5240  
    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 
    … …  
    2528   USE wrk_nemo       ! Memory allocation 
    2629   USE timing         ! Timing 
     30   USE iom            ! I/O manager 
    2731 
    2832   IMPLICIT NONE 
    … …  
    7377      NAMELIST/namc1d_dyndmp/ ln_dyndmp 
    7478      INTEGER :: ios 
     79      INTEGER :: imask 
    7580      !!---------------------------------------------------------------------- 
    7681 
    … …  
    9196         WRITE(numout,*) '      add a damping term or not       ln_dyndmp = ', ln_dyndmp 
    9297         WRITE(numout,*) '   Namelist namtra_dmp    : Set damping parameters' 
    93          WRITE(numout,*) '      horizontal damping option       nn_hdmp   = ', nn_hdmp 
    94          WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp, '(non-C1D zoom: forced to 0)' 
    95          WRITE(numout,*) '      surface time scale (days)       rn_surf   = ', rn_surf 
    96          WRITE(numout,*) '      bottom time scale (days)        rn_bot    = ', rn_bot 
    97          WRITE(numout,*) '      depth of transition (meters)    rn_dep    = ', rn_dep 
    98          WRITE(numout,*) '      create a damping.coeff file     nn_file   = ', nn_file 
     98         WRITE(numout,*) '      Apply relaxation   or not       ln_tradmp = ', ln_tradmp 
     99         WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp 
     100         WRITE(numout,*) '      Damping file name               cn_resto  = ', cn_resto 
    99101         WRITE(numout,*) 
    100102      ENDIF 
    … …  
    104106         IF( dyn_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dyn_dmp_init: unable to allocate arrays' ) 
    105107         ! 
    106 #if ! defined key_c1d 
    107          SELECT CASE ( nn_hdmp )             !==   control print of horizontal option   ==! 
    108          CASE (  -1  )   ;   IF(lwp) WRITE(numout,*) '   momentum damping in the Med & Red seas only' 
    109          CASE ( 1:90 )   ;   IF(lwp) WRITE(numout,*) '   momentum damping poleward of', nn_hdmp, ' degrees' 
    110          CASE DEFAULT 
    111             WRITE(ctmp1,*) '          bad flag value for nn_hdmp = ', nn_hdmp 
    112             CALL ctl_stop(ctmp1) 
    113          END SELECT 
    114          ! 
    115 #endif 
    116108         SELECT CASE ( nn_zdmp )             !==   control print of vertical option   ==! 
    117109         CASE ( 0    )   ;   IF(lwp) WRITE(numout,*) '   momentum damping throughout the water column' 
    … …  
    130122         utrdmp(:,:,:) = 0._wp               ! internal damping trends 
    131123         vtrdmp(:,:,:) = 0._wp 
    132          !                                   !==   Damping coefficients calculation:                           ==! 
    133          !                                   !==   use tradmp.F90 subroutines dtacof, dtacof_zoom and cofdis   ==! 
    134          !                                   !!!   NOTE: these need to be altered for use in this module if  
    135          !                                   !!!       they are to be used outside the C1D context  
    136          !                                   !!!       (use of U,V grid variables) 
    137          IF( lzoom .AND. .NOT. lk_c1d ) THEN   ;   CALL dtacof_zoom( resto_uv ) 
    138          ELSE               ;   CALL dtacof( nn_hdmp, rn_surf, rn_bot, rn_dep, nn_file, 'DYN', resto_uv ) 
    139          ENDIF 
    140          ! 
     124         ! 
     125         !Read in mask from file 
     126         CALL iom_open ( cn_resto, imask) 
     127         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto) 
     128         CALL iom_close( imask ) 
    141129      ENDIF 
    142130      ! 

  • branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

    r5239 r5240  
    8282      CALL wrk_alloc( jpi , jpj , jpk        , zrhd      , zrhop    ) 
    8383      CALL wrk_alloc( jpi , jpj , jpk , jpts , ztsn                 ) 
    84  
    85       CALL iom_put( 'cellthc', fse3t(:,:,:) ) 
    8684 
    8785      zarea_ssh(:,:) = area(:,:) * sshn(:,:) 

  • branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r5239 r5240  
    143143      ENDIF 
    144144 
    145       IF( lk_vvl ) THEN 
    146          z3d(:,:,:) = tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) 
    147          CALL iom_put( "toce" , z3d                        )   ! heat content 
     145      IF( .NOT.lk_vvl ) THEN 
     146         CALL iom_put( "e3t" , fse3t_n(:,:,:) ) 
     147         CALL iom_put( "e3u" , fse3u_n(:,:,:) ) 
     148         CALL iom_put( "e3v" , fse3v_n(:,:,:) ) 
     149         CALL iom_put( "e3w" , fse3w_n(:,:,:) ) 
     150      ENDIF 
     151       
     152      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature 
     153      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature 
     154      IF ( iom_use("sbt") ) THEN 
    148155         DO jj = 1, jpj 
    149156            DO ji = 1, jpi 
    150                z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * fse3t_n(ji,jj,mikt(ji,jj)) 
    151             END DO 
    152          END DO   
    153          CALL iom_put( "sst"  , z2d(:,:)                 )   ! sea surface heat content       
     157               z2d(ji,jj) = tsn(ji,jj,MAX(mbathy(ji,jj),1),jp_tem) 
     158            END DO 
     159         END DO 
     160         CALL iom_put( "sbt", z2d )                ! bottom temperature 
     161      ENDIF 
     162       
     163      CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity 
     164      CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity 
     165      IF ( iom_use("sbs") ) THEN 
    154166         DO jj = 1, jpj 
    155167            DO ji = 1, jpi 
    156                z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem)**2 * fse3t_n(ji,jj,mikt(ji,jj)) 
    157             END DO 
    158          END DO   
    159          CALL iom_put( "sst2" , z2d(:,:)      )   ! sea surface content of squared temperature 
    160          z3d(:,:,:) = tsn(:,:,:,jp_sal) * fse3t_n(:,:,:)             
    161          CALL iom_put( "soce" , z3d                        )   ! salinity content 
     168               z2d(ji,jj) = tsn(ji,jj,MAX(mbathy(ji,jj),1),jp_sal) 
     169            END DO 
     170         END DO 
     171         CALL iom_put( "sbs", z2d )                ! bottom salinity 
     172      ENDIF 
     173          
     174      CALL iom_put( "uoce", un(:,:,:)         )    ! 3D i-current 
     175      CALL iom_put(  "ssu", un(:,:,1)         )    ! surface i-current 
     176      IF ( iom_use("sbu") ) THEN 
    162177         DO jj = 1, jpj 
    163178            DO ji = 1, jpi 
    164                z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * fse3t_n(ji,jj,mikt(ji,jj)) 
    165             END DO 
    166          END DO   
    167          CALL iom_put( "sss"  , z2d(:,:)                 )   ! sea surface salinity content 
     179               z2d(ji,jj) = un(ji,jj,MAX(mbathy(ji,jj),1)) 
     180            END DO 
     181         END DO 
     182         CALL iom_put( "sbu", z2d )                ! bottom i-current 
     183      ENDIF 
     184       
     185      CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current 
     186      CALL iom_put(  "ssv", vn(:,:,1)         )    ! surface j-current 
     187      IF ( iom_use("sbv") ) THEN 
    168188         DO jj = 1, jpj 
    169189            DO ji = 1, jpi 
    170                z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal)**2 * fse3t_n(ji,jj,mikt(ji,jj)) 
    171             END DO 
    172          END DO   
    173          CALL iom_put( "sss2" , z2d(:,:)                 )   ! sea surface content of squared salinity 
    174       ELSE 
    175          CALL iom_put( "toce" , tsn(:,:,:,jp_tem)        )   ! temperature 
    176          IF ( iom_use("sst") ) THEN 
    177             DO jj = 1, jpj 
    178                DO ji = 1, jpi 
    179                   z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) 
    180                END DO 
    181             END DO 
    182             CALL iom_put( "sst"  , z2d(:,:)            ) ! sea surface temperature 
    183          ENDIF 
    184          IF ( iom_use("sst2") )   CALL iom_put( "sst2" , z2d(:,:) * z2d(:,:) ) ! square of sea surface temperature 
    185          CALL iom_put( "soce" , tsn(:,:,:,jp_sal)          )   ! salinity 
    186          IF ( iom_use("sss") ) THEN 
    187             DO jj = 1, jpj 
    188                DO ji = 1, jpi 
    189                   z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) 
    190                END DO 
    191             END DO 
    192             CALL iom_put( "sss"  , z2d(:,:)            ) ! sea surface salinity 
    193          ENDIF 
    194          CALL iom_put( "sss2" , z2d(:,:) * z2d(:,:) ) ! square of sea surface salinity 
    195       END IF 
    196       IF( lk_vvl .AND. (.NOT. ln_dynadv_vec) ) THEN 
    197          CALL iom_put( "uoce" , umask(:,:,:) * un(:,:,:) * fse3u_n(:,:,:) )    ! i-transport 
    198          CALL iom_put( "voce" , vmask(:,:,:) * vn(:,:,:) * fse3v_n(:,:,:) )    ! j-transport 
    199       ELSE 
    200          CALL iom_put( "uoce" , umask(:,:,:) * un(:,:,:)                  )    ! i-current 
    201          CALL iom_put( "voce" , vmask(:,:,:) * vn(:,:,:)                  )    ! j-current 
    202          IF ( iom_use("ssu") ) THEN 
    203             DO jj = 1, jpj 
    204                DO ji = 1, jpi 
    205                   z2d(ji,jj) = un(ji,jj,miku(ji,jj)) 
    206                END DO 
    207             END DO 
    208             CALL iom_put( "ssu"   , z2d                                    )    ! i-current 
    209          ENDIF 
    210          IF ( iom_use("ssv") ) THEN 
    211             DO jj = 1, jpj 
    212                DO ji = 1, jpi 
    213                   z2d(ji,jj) = vn(ji,jj,mikv(ji,jj)) 
    214                END DO 
    215             END DO 
    216             CALL iom_put( "ssv"   , z2d                                    )    ! j-current 
    217          ENDIF 
    218       ENDIF 
    219       CALL iom_put(    "avt"  , avt                        )    ! T vert. eddy diff. coef. 
    220       CALL iom_put(    "avm"  , avmu                       )    ! T vert. eddy visc. coef. 
    221       IF( lk_zdftke ) THEN   
    222          CALL iom_put( "tke"      , en                               )    ! TKE budget: Turbulent Kinetic Energy   
    223          CALL iom_put( "tke_niw"  , e_niw                            )    ! TKE budget: Near-inertial waves   
    224       ENDIF 
    225       IF( lk_zdfddm ) THEN 
    226          CALL iom_put( "avs" , fsavs(:,:,:)                          )    ! S vert. eddy diff. coef. 
    227       ENDIF 
    228  
    229       IF ( iom_use("sstgrad2") .OR. iom_use("sstgrad2") ) THEN 
     190               z2d(ji,jj) = vn(ji,jj,MAX(mbathy(ji,jj),1)) 
     191            END DO 
     192         END DO 
     193         CALL iom_put( "sbv", z2d )                ! bottom j-current 
     194      ENDIF 
     195 
     196      CALL iom_put( "avt" , avt                        )    ! T vert. eddy diff. coef. 
     197      CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef. 
     198      IF( lk_zdftke ) THEN    
     199         CALL iom_put( "tke"      , en                               )    ! TKE budget: Turbulent Kinetic Energy    
     200         CALL iom_put( "tke_niw"  , e_niw                            )    ! TKE budget: Near-inertial waves    
     201      ENDIF  
     202      CALL iom_put( "avs" , fsavs(:,:,:)               )    ! S vert. eddy diff. coef. (useful only with key_zdfddm) 
     203 
     204      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN 
    230205         DO jj = 2, jpjm1                                    ! sst gradient 
    231206            DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    239214         CALL lbc_lnk( z2d, 'T', 1. ) 
    240215         CALL iom_put( "sstgrad2",  z2d               )    ! square of module of sst gradient 
    241          !CDIR NOVERRCHK< 
    242216         z2d(:,:) = SQRT( z2d(:,:) ) 
    243217         CALL iom_put( "sstgrad" ,  z2d               )    ! module of sst gradient 
    … …  
    248222         z2d(:,:)  = 0._wp  
    249223         DO jk = 1, jpkm1 
    250             DO jj = 2, jpjm1 
    251                DO ji = fs_2, fs_jpim1   ! vector opt. 
     224            DO jj = 1, jpj 
     225               DO ji = 1, jpi 
    252226                  z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 
    253227               END DO 
    254228            END DO 
    255229         END DO 
    256          CALL lbc_lnk( z2d, 'T', 1. ) 
    257230         CALL iom_put( "heatc", (rau0 * rcp) * z2d )    ! vertically integrated heat content (J/m2) 
    258231      ENDIF 
    … …  
    261234         z2d(:,:)  = 0._wp  
    262235         DO jk = 1, jpkm1 
    263             DO jj = 2, jpjm1 
    264                DO ji = fs_2, fs_jpim1   ! vector opt. 
     236            DO jj = 1, jpj 
     237               DO ji = 1, jpi 
    265238                  z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk) 
    266239               END DO 
    267240            END DO 
    268241         END DO 
    269          CALL lbc_lnk( z2d, 'T', 1. ) 
    270242         CALL iom_put( "saltc", rau0 * z2d )   ! vertically integrated salt content (PSU*kg/m2) 
    271243      ENDIF 

  • branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r5239 r5240  
    588588      INTEGER, INTENT( in )               :: kt       ! time step 
    589589      !! * Local declarations 
    590       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def 
    591590      INTEGER                             :: ji,jj,jk       ! dummy loop indices 
    592591      !!---------------------------------------------------------------------- 
    593592 
    594593      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp') 
    595       ! 
    596       CALL wrk_alloc( jpi, jpj, jpk, z_e3t_def                ) 
    597594      ! 
    598595      IF( kt == nit000 )   THEN 
    … …  
    679676      ! Write outputs 
    680677      ! ============= 
    681       z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 
    682       CALL iom_put( "cellthc" , fse3t_n  (:,:,:) ) 
     678      CALL iom_put(     "e3t" , fse3t_n  (:,:,:) ) 
     679      CALL iom_put(     "e3u" , fse3u_n  (:,:,:) ) 
     680      CALL iom_put(     "e3v" , fse3v_n  (:,:,:) ) 
     681      CALL iom_put(     "e3w" , fse3w_n  (:,:,:) ) 
    683682      CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) ) 
    684       CALL iom_put( "e3tdef"  , z_e3t_def(:,:,:) ) 
     683      IF( iom_use("e3tdef") )   & 
     684         CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
    685685 
    686686      ! write restart file 
    687687      ! ================== 
    688688      IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' ) 
    689       ! 
    690       CALL wrk_dealloc( jpi, jpj, jpk, z_e3t_def ) 
    691689      ! 
    692690      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_swp') 

  • branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5239 r5240  
    365365      INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices 
    366366      INTEGER  ::   inum                      ! temporary logical unit 
     367      INTEGER  ::   ierror                    ! error flag 
    367368      INTEGER  ::   ii_bump, ij_bump, ih      ! bump center position 
    368369      INTEGER  ::   ii0, ii1, ij0, ij1, ik    ! local indices 
    369370      REAL(wp) ::   r_bump , h_bump , h_oce   ! bump characteristics  
    370371      REAL(wp) ::   zi, zj, zh, zhmin         ! local scalars 
    371       INTEGER , POINTER, DIMENSION(:,:) ::   idta   ! global domain integer data 
    372       REAL(wp), POINTER, DIMENSION(:,:) ::   zdta   ! global domain scalar data 
     372      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   idta   ! global domain integer data 
     373      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdta   ! global domain scalar data 
    373374      !!---------------------------------------------------------------------- 
    374375      ! 
    375376      IF( nn_timing == 1 )  CALL timing_start('zgr_bat') 
    376       ! 
    377       CALL wrk_alloc( jpidta, jpjdta, idta ) 
    378       CALL wrk_alloc( jpidta, jpjdta, zdta ) 
    379377      ! 
    380378      IF(lwp) WRITE(numout,*) 
    … …  
    385383         !                                            ! ================== ! 
    386384         !                                            ! global domain level and meter bathymetry (idta,zdta) 
     385         ! 
     386         ALLOCATE( idta(jpidta,jpjdta), STAT=ierror ) 
     387         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' ) 
     388         ALLOCATE( zdta(jpidta,jpjdta), STAT=ierror ) 
     389         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' ) 
    387390         ! 
    388391         IF( ntopo == 0 ) THEN                        ! flat basin 
    … …  
    489492            WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
    490493         END IF 
     494         ! 
     495         DEALLOCATE( idta, zdta ) 
    491496         ! 
    492497         !                                            ! ================ ! 
    … …  
    593598      ENDIF 
    594599      ! 
    595       CALL wrk_dealloc( jpidta, jpjdta, idta ) 
    596       CALL wrk_dealloc( jpidta, jpjdta, zdta ) 
    597       ! 
    598600      IF( nn_timing == 1 )  CALL timing_stop('zgr_bat') 
    599601      ! 

  • branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90

    r5239 r5240  
    116116         DO jj = 2, jpjm1                          ! laplacian 
    117117            DO ji = fs_2, fs_jpim1   ! vector opt. 
    118                zlu_uu(ji,jj,jk,1) = ( ub (ji+1,jj,jk)-2.*ub (ji,jj,jk)+ub (ji-1,jj,jk) ) * umask(ji,jj,jk) 
    119                zlv_vv(ji,jj,jk,1) = ( vb (ji,jj+1,jk)-2.*vb (ji,jj,jk)+vb (ji,jj-1,jk) ) * vmask(ji,jj,jk)  
    120                zlu_uv(ji,jj,jk,1) = ( ub (ji,jj+1,jk)-2.*ub (ji,jj,jk)+ub (ji,jj-1,jk) ) * umask(ji,jj,jk) 
    121                zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj,jk)-2.*vb (ji,jj,jk)+vb (ji-1,jj,jk) ) * vmask(ji,jj,jk) 
    122                ! 
    123                zlu_uu(ji,jj,jk,2) = ( zfu(ji+1,jj,jk)-2.*zfu(ji,jj,jk)+zfu(ji-1,jj,jk) ) * umask(ji,jj,jk) 
    124                zlv_vv(ji,jj,jk,2) = ( zfv(ji,jj+1,jk)-2.*zfv(ji,jj,jk)+zfv(ji,jj-1,jk) ) * vmask(ji,jj,jk) 
    125                zlu_uv(ji,jj,jk,2) = ( zfu(ji,jj+1,jk)-2.*zfu(ji,jj,jk)+zfu(ji,jj-1,jk) ) * umask(ji,jj,jk) 
    126                zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj,jk)-2.*zfv(ji,jj,jk)+zfv(ji-1,jj,jk) ) * vmask(ji,jj,jk) 
    127             END DO 
    128          END DO 
    129       END DO 
    130 !!gm BUG !!!  just below this should be +1 in all the communications 
    131 !      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', -1.)   ;   CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', -1.) 
    132 !      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', -1.)   ;   CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', -1.) 
    133 !      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', -1.)   ;   CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', -1.) 
    134 !      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', -1.)   ;   CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', -1.) 
    135 ! 
    136 !!gm corrected: 
     118               ! 
     119               zlu_uu(ji,jj,jk,1) = ( ub (ji+1,jj  ,jk) - 2.*ub (ji,jj,jk) + ub (ji-1,jj  ,jk) ) * umask(ji,jj,jk) 
     120               zlv_vv(ji,jj,jk,1) = ( vb (ji  ,jj+1,jk) - 2.*vb (ji,jj,jk) + vb (ji  ,jj-1,jk) ) * vmask(ji,jj,jk) 
     121               zlu_uv(ji,jj,jk,1) = ( ub (ji  ,jj+1,jk) - ub (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   & 
     122                  &               - ( ub (ji  ,jj  ,jk) - ub (ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk) 
     123               zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj  ,jk) - vb (ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   & 
     124                  &               - ( vb (ji  ,jj  ,jk) - vb (ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk) 
     125               ! 
     126               zlu_uu(ji,jj,jk,2) = ( zfu(ji+1,jj  ,jk) - 2.*zfu(ji,jj,jk) + zfu(ji-1,jj  ,jk) ) * umask(ji,jj,jk) 
     127               zlv_vv(ji,jj,jk,2) = ( zfv(ji  ,jj+1,jk) - 2.*zfv(ji,jj,jk) + zfv(ji  ,jj-1,jk) ) * vmask(ji,jj,jk) 
     128               zlu_uv(ji,jj,jk,2) = ( zfu(ji  ,jj+1,jk) - zfu(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   & 
     129                  &               - ( zfu(ji  ,jj  ,jk) - zfu(ji  ,jj-1,jk) ) * fmask(ji  ,jj-1,jk) 
     130               zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj  ,jk) - zfv(ji  ,jj  ,jk) ) * fmask(ji  ,jj  ,jk)   & 
     131                  &               - ( zfv(ji  ,jj  ,jk) - zfv(ji-1,jj  ,jk) ) * fmask(ji-1,jj  ,jk) 
     132            END DO 
     133         END DO 
     134      END DO 
    137135      CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', 1. ) 
    138136      CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', 1. )   ;   CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', 1. ) 
    139137      CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', 1. ) 
    140138      CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', 1. )   ;   CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', 1. )  
    141 !!gm end 
    142139       
    143140      !                                      ! ====================== ! 

  • branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r5239 r5240  
    9797      ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) 
    9898 
    99       IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &  
    100                              &      ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 
     99      IF( ln_dynvor_een .or. ln_dynvor_een_old ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &  
     100                                                    &      ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 
    101101 
    102102      dyn_spg_ts_alloc = MAXVAL(ierr(:)) 
    … …  
    218218      ! 
    219219      IF ( kt == nit000 .OR. lk_vvl ) THEN 
    220          IF ( ln_dynvor_een ) THEN 
     220         IF ( ln_dynvor_een_old ) THEN 
     221            DO jj = 1, jpjm1 
     222               DO ji = 1, jpim1 
     223                  zwz(ji,jj) =   ( ht(ji  ,jj+1) + ht(ji+1,jj+1) +                    & 
     224                        &          ht(ji  ,jj  ) + ht(ji+1,jj  )   ) / 4._wp   
     225                  IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = 1._wp / zwz(ji,jj) 
     226               END DO 
     227            END DO 
     228            CALL lbc_lnk( zwz, 'F', 1._wp ) 
     229            zwz(:,:) = ff(:,:) * zwz(:,:) 
     230 
     231            ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 
     232            DO jj = 2, jpj 
     233               DO ji = fs_2, jpi   ! vector opt. 
     234                  ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1) 
     235                  ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  ) 
     236                  ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1) 
     237                  ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  ) 
     238               END DO 
     239            END DO 
     240         ELSE IF ( ln_dynvor_een ) THEN 
    221241            DO jj = 1, jpjm1 
    222242               DO ji = 1, jpim1 
    … …  
    339359         END DO 
    340360         ! 
    341       ELSEIF ( ln_dynvor_een ) THEN                    ! enstrophy and energy conserving scheme 
     361      ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN  ! enstrophy and energy conserving scheme 
    342362         DO jj = 2, jpjm1 
    343363            DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    687707            END DO 
    688708            ! 
    689          ELSEIF ( ln_dynvor_een ) THEN                    !==  energy and enstrophy conserving scheme  ==! 
     709         ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN !==  energy and enstrophy conserving scheme  ==! 
    690710            DO jj = 2, jpjm1 
    691711               DO ji = fs_2, fs_jpim1   ! vector opt. 

  • branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90

    r5239 r5240  
    5151   LOGICAL, PUBLIC ::   ln_dynvor_mix   !: mixed scheme 
    5252   LOGICAL, PUBLIC ::   ln_dynvor_een   !: energy and enstrophy conserving scheme 
     53   LOGICAL, PUBLIC ::   ln_dynvor_een_old !: energy and enstrophy conserving scheme (original formulation) 
    5354 
    5455   INTEGER ::   nvor = 0   ! type of vorticity trend used 
    … …  
    596597 
    597598      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points) 
    598          DO jk = 1, jpk 
    599             DO jj = 1, jpjm1 
    600                DO ji = 1, jpim1 
    601                   ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
    602                      &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) 
    603                   zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   & 
    604                      &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) ) 
    605                   IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = zmsk / ze3 
    606                END DO 
    607             END DO 
    608          END DO 
     599 
     600         IF( ln_dynvor_een_old ) THEN ! original formulation 
     601            DO jk = 1, jpk 
     602               DO jj = 1, jpjm1 
     603                  DO ji = 1, jpim1 
     604                     ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
     605                        &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) 
     606                     IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = 4.0_wp / ze3 
     607                  END DO 
     608               END DO 
     609            END DO 
     610         ELSE ! new formulation from NEMO 3.6 
     611            DO jk = 1, jpk 
     612               DO jj = 1, jpjm1 
     613                  DO ji = 1, jpim1 
     614                     ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   & 
     615                        &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) 
     616                     zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   & 
     617                        &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) ) 
     618                     IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = zmsk / ze3 
     619                  END DO 
     620               END DO 
     621            END DO 
     622         ENDIF 
     623 
    609624         CALL lbc_lnk( ze3f, 'F', 1. ) 
    610625      ENDIF 
    … …  
    705720      INTEGER ::   ios             ! Local integer output status for namelist read 
    706721      !! 
    707       NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een 
     722      NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, ln_dynvor_een_old 
    708723      !!---------------------------------------------------------------------- 
    709724 
    … …  
    726741         WRITE(numout,*) '           mixed enstrophy/energy conserving scheme   ln_dynvor_mix = ', ln_dynvor_mix 
    727742         WRITE(numout,*) '           enstrophy and energy conserving scheme     ln_dynvor_een = ', ln_dynvor_een 
     743         WRITE(numout,*) '           enstrophy and energy conserving scheme (old) ln_dynvor_een_old= ', ln_dynvor_een_old 
    728744      ENDIF 
    729745 
    … …  
    749765      IF( ln_dynvor_mix )   ioptio = ioptio + 1 
    750766      IF( ln_dynvor_een )   ioptio = ioptio + 1 
     767      IF( ln_dynvor_een_old )   ioptio = ioptio + 1 
    751768      IF( lk_esopa      )   ioptio =          1 
    752769 
    … …  
    757774      IF( ln_dynvor_ens )   nvor =  1 
    758775      IF( ln_dynvor_mix )   nvor =  2 
    759       IF( ln_dynvor_een )   nvor =  3 
     776      IF( ln_dynvor_een .or. ln_dynvor_een_old )   nvor =  3 
    760777      IF( lk_esopa      )   nvor = -1 
    761778       

  • branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90

    r5239 r5240  
    12021202         CASE('T')   ;   zmask(:,:,:)       = tmask(:,:,:) 
    12031203         CASE('U')   ;   zmask(2:jpim1,:,:) = tmask(2:jpim1,:,:) + tmask(3:jpi,:,:)   ;   CALL lbc_lnk( zmask, 'U', 1. ) 
    1204          CASE('V')   ;   zmask(:,2:jpjm1,:) = tmask(:,2:jpjm1,:) + tmask(:,3:jpi,:)   ;   CALL lbc_lnk( zmask, 'V', 1. ) 
     1204         CASE('V')   ;   zmask(:,2:jpjm1,:) = tmask(:,2:jpjm1,:) + tmask(:,3:jpj,:)   ;   CALL lbc_lnk( zmask, 'V', 1. ) 
    12051205         CASE('W')   ;   zmask(:,:,2:jpk  ) = tmask(:,:,1:jpkm1) + tmask(:,:,2:jpk)   ;   zmask(:,:,1) = tmask(:,:,1) 
    12061206         END SELECT 

  • branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/IOM/prtctl.F90

    r5239 r5240  
    164164         ENDIF 
    165165 
    166          IF ( clinfo3 == 'tra' )  THEN 
    167              zvctl1 = t_ctll(jn) 
    168              zvctl2 = s_ctll(jn) 
    169          ELSEIF ( clinfo3 == 'dyn' )   THEN 
    170              zvctl1 = u_ctll(jn) 
    171              zvctl2 = v_ctll(jn) 
     166         IF( PRESENT(clinfo3)) THEN 
     167            IF ( clinfo3 == 'tra' )  THEN 
     168               zvctl1 = t_ctll(jn) 
     169               zvctl2 = s_ctll(jn) 
     170            ELSEIF ( clinfo3 == 'dyn' )   THEN 
     171               zvctl1 = u_ctll(jn) 
     172               zvctl2 = v_ctll(jn) 
     173            ENDIF 
    172174         ENDIF 
    173175 

  • branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn_c2d.h90

    r5239 r5240  
    146146      INTEGER  ::   inum, iim, ijm            ! local integers 
    147147      INTEGER  ::   ifreq, il1, il2, ij, ii 
    148       INTEGER  ::   ijpt0,ijpt1 
     148      INTEGER  ::   ijpt0,ijpt1, ierror 
    149149      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk 
    150150      CHARACTER (len=15) ::   clexp 
    151       INTEGER, POINTER, DIMENSION(:,:)  :: icof 
    152       INTEGER, POINTER, DIMENSION(:,:)  :: idata 
     151      INTEGER,     POINTER, DIMENSION(:,:)  :: icof 
     152      INTEGER, ALLOCATABLE, DIMENSION(:,:)  :: idata 
    153153      !!---------------------------------------------------------------------- 
    154154      !                                 
    155155      CALL wrk_alloc( jpi   , jpj   , icof  ) 
    156       CALL wrk_alloc( jpidta, jpjdta, idata ) 
    157156      ! 
    158157      IF(lwp) WRITE(numout,*) 
    … …  
    234233         ! ===================== equatorial strip (20N-20S) defined at t-points 
    235234          
     235         ALLOCATE( idata(jpidta,jpjdta), STAT=ierror ) 
     236         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_c2d_orca: unable to allocate idata array' ) 
     237         ! 
    236238         CALL ctl_opn( inum, 'ahmcoef', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
    237239         READ(inum,9101) clexp, iim, ijm 
    … …  
    2692719201     FORMAT(3x,13(i3,12x)) 
    2702729202     FORMAT(i3,41i3) 
    271  
     273          
     274         DEALLOCATE(idata) 
    272275 
    273276         ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator) 
    … …  
    346349      ! 
    347350      CALL wrk_dealloc( jpi   , jpj   , icof  ) 
    348       CALL wrk_dealloc( jpidta, jpjdta, idata ) 
    349351      ! 
    350352   END SUBROUTINE ldf_dyn_c2d_orca 
    … …  
    374376      INTEGER ::   iim, ijm 
    375377      INTEGER ::   ifreq, il1, il2, ij, ii 
    376       INTEGER ::   ijpt0,ijpt1 
     378      INTEGER ::   ijpt0,ijpt1, ierror 
    377379      REAL(wp) ::   zahmeq, zcoft, zcoff, zmsk, zam20s 
    378380      CHARACTER (len=15) ::   clexp 
    379       INTEGER, POINTER, DIMENSION(:,:)  :: icof 
    380       INTEGER, POINTER, DIMENSION(:,:)  :: idata 
     381      INTEGER,     POINTER, DIMENSION(:,:)  :: icof 
     382      INTEGER, ALLOCATABLE, DIMENSION(:,:)  :: idata 
    381383      !!---------------------------------------------------------------------- 
    382384      !                                 
    383385      CALL wrk_alloc( jpi   , jpj   , icof  ) 
    384       CALL wrk_alloc( jpidta, jpjdta, idata ) 
    385386      !                                 
    386387 
    … …  
    464465         ! ===================== equatorial strip (20N-20S) defined at t-points 
    465466          
     467         ALLOCATE( idata(jpidta,jpjdta), STAT=ierror ) 
     468         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_c2d_orca_R1: unable to allocate idata array' ) 
     469         ! 
    466470         CALL ctl_opn( inum, 'ahmcoef', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',   & 
    467471            &           1, numout, lwp ) 
    … …  
    5015059201     FORMAT(3x,13(i3,12x)) 
    5025069202     FORMAT(i3,41i3) 
    503  
     507          
     508         DEALLOCATE(idata) 
    504509 
    505510         ! Set ahm1 and ahm2  ( T- and F- points) (used for laplacian operator) 
    … …  
    583588      ! 
    584589      CALL wrk_dealloc( jpi   , jpj   , icof  ) 
    585       CALL wrk_dealloc( jpidta, jpjdta, idata ) 
    586590      ! 
    587591   END SUBROUTINE ldf_dyn_c2d_orca_R1 

  • branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r5239 r5240  
    848848      rgt33 = 0.5_wp + SIGN( 0.5_wp, (zw10 - 33._wp) )   ! If zw10 < 33. => 0, else => 1   
    849849      cd_neutral_10m = 1.e-3 * ( & 
    850          &       (rgt33 + 1._wp)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33. 
     850         &       (1._wp - rgt33)*( 2.7_wp/zw10 + 0.142_wp + zw10/13.09_wp - 3.14807E-10*zw10**6) & ! zw10< 33. 
    851851         &      + rgt33         *      2.34   )                                                    ! zw10 >= 33. 
    852852      ! 

  • branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90

    r5239 r5240  
    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.