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 9149 – NEMO

Changeset 9149


Ignore:
Timestamp:
2017-12-20T16:38:27+01:00 (6 years ago)
Author:
jchanut
Message:

NESTING TOOLS
Add new e3t definition, #1981
Correct transition weights, #1998
Correct arithmetic and median average interpolations, #1999
Changed number of coarse grid bathymetry points inside child grid domain (2 now, instead of 3)

Location:
branches/2017/dev_merge_2017/NEMOGCM/TOOLS/NESTING
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_merge_2017/NEMOGCM/TOOLS/NESTING/agulhas

    r6819 r9149  
    3232 
    3333&vertical_grid     
     34    ln_e3_dep = .true. 
    3435    ppkth = 21.4333619793800 
    3536    ppacr = 3 
  • branches/2017/dev_merge_2017/NEMOGCM/TOOLS/NESTING/src/agrif_connect_topo.f90

    r5656 r9149  
    188188    ! 
    189189    gdepw(1) = 0.0 
     190    ! 
     191    IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]    
     192       ! 
     193       DO i = 1, N-1  
     194          e3t(i) = gdepw(i+1)-gdepw(i) 
     195       END DO 
     196       e3t(N) = e3t(N-1) 
     197 
     198       DO i = 2, N  
     199          e3w(i) = gdept(i) - gdept(i-1) 
     200       END DO 
     201       e3w(1  ) = 2. * (gdept(1) - gdepw(1)) 
     202    END IF   
     203    ! 
    190204    zmax = gdepw(N) + e3t(N) 
    191205    IF( rn_hmin < 0. ) THEN  ;   i = - INT( rn_hmin )                                  ! from a nb of level 
     
    342356    ! 
    343357    gdepw(1) = 0.0   
     358    ! 
     359    IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]    
     360       ! 
     361       DO i = 1, N-1 
     362          e3t(i) = gdepw(i+1)-gdepw(i) 
     363       END DO 
     364       e3t(N) = e3t(N-1) 
     365 
     366       DO i = 2, N 
     367          e3w(i) = gdept(i) - gdept(i-1) 
     368       END DO 
     369       e3w(1  ) = 2. * (gdept(1) - gdepw(1)) 
     370    END IF 
    344371    ! 
    345372    IF(.NOT. ASSOCIATED(Grid%bathy_meter)) THEN 
     
    838865    ENDIF 
    839866    ! 
     867    gdepw(1)=0.  
     868    IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]    
     869       ! 
     870       DO i = 1, N-1 
     871          e3t(i) = gdepw(i+1)-gdepw(i) 
     872       END DO 
     873       e3t(N) = e3t(N-1) 
     874 
     875       DO i = 2, N 
     876          e3w(i) = gdept(i) - gdept(i-1) 
     877       END DO 
     878       e3w(1  ) = 2. * (gdept(1) - gdepw(1)) 
     879    END IF 
     880    ! 
    840881    IF( rn_hmin < 0. ) THEN  ;   i = - INT( rn_hmin )                                  ! from a nb of level 
    841882    ELSE                     ;   i = MINLOC( gdepw, mask = gdepw > rn_hmin, dim = 1 )  ! from a depth 
  • branches/2017/dev_merge_2017/NEMOGCM/TOOLS/NESTING/src/agrif_create_bathy.f90

    r8641 r9149  
    5353  REAL*8, DIMENSION(:,:),POINTER :: gdepw_ps_interp => NULL()  
    5454  REAL*8, DIMENSION(:,:),POINTER :: save_gdepw,rx,ry,maskedtopo 
     55  REAL*8, DIMENSION(:,:),POINTER :: wtab  => NULL() 
    5556  REAL*8  :: Cell_lonmin,Cell_lonmax,Cell_latmin,Cell_latmax,wghts 
    5657  LOGICAL :: Pacifique=.FALSE. 
     
    239240        status = Read_Bathymeter(TRIM(parent_bathy_meter),G0) 
    240241        IF(Pacifique) THEN 
    241            WHERE(G0%nav_lon < 0.001)  
     242           WHERE(G0%nav_lon < 0.0001)  
    242243              G0%nav_lon = G0%nav_lon + 360. 
    243244           END WHERE 
     
    302303              DO WHILE( G0%nav_lon(iimax,1) <= Cell_lonmax )  
    303304                 iimax = iimax + 1 
     305              iimax = MIN( iimax,SIZE(G0%bathy_meter,1)) 
    304306              ENDDO 
    305307              !                                
     
    307309              DO WHILE( G0%nav_lat(iimax,jjmax) <= Cell_latmax )  
    308310                 jjmax = jjmax + 1 
     311              jjmax = MIN( jjmax,SIZE(G0%bathy_meter,2)) 
     312 
    309313              ENDDO 
    310314              ! 
     
    439443        status = Read_coordinates(TRIM(parent_coordinate_file),G0) 
    440444        !------------------------                   
    441  
    442445        IF (.NOT.ASSOCIATED(G0%gdepw_ps)) & 
    443446             ALLOCATE(G0%gdepw_ps(SIZE(G0%bathy_meter,1),SIZE(G0%bathy_meter,2))) 
     
    468471        WRITE(*,*) ' linear connection on ',nb_connection_pts,'coarse grid points' 
    469472 
    470         connectionsize = 3 + nb_connection_pts  
     473        connectionsize = 2 + nb_connection_pts  
    471474        !             
    472475        gdepw_ps_interp = 0. 
     
    482485        ! 
    483486        ! 
     487        IF (.NOT.ASSOCIATED(wtab)) & 
     488             ALLOCATE(wtab(SIZE(G1%bathy_meter,1),SIZE(G1%bathy_meter,2))) 
     489        wtab(:,:) = 0. 
    484490        wghts = 1. 
    485491        DO ji = boundary + 1 , boundary + nb_connection_pts * irafx 
    486            G1%gdepw_ps(ji,boundary+1:nyfin-boundary) =                                          & 
    487                 (1.-wghts) * G1%gdepw_ps(ji,boundary+1:nyfin-boundary) +                          & 
    488                 wghts * gdepw_ps_interp(ji,boundary+1:nyfin-boundary) 
    489            wghts = wghts - (1. / (nb_connection_pts*irafx - 1. ) )  
    490         ENDDO 
    491  
     492            wghts = wghts - (1. / (nb_connection_pts*irafx - 1. ) ) 
     493            DO jj=boundary+1,nyfin-boundary 
     494               wtab(ji,jj) = MAX(wghts, wtab(ji,jj))   
     495            END DO 
     496        END DO  
     497       
    492498        wghts = 1. 
    493         DO ji = nxfin - boundary , nxfin - boundary - nb_connection_pts * irafx + 1 ,-1  
    494            G1%gdepw_ps(ji,boundary+1:nyfin-boundary) =                                            & 
    495                 (1. - wghts) * G1%gdepw_ps(ji,boundary+1:nyfin-boundary) +                          & 
    496                 wghts * gdepw_ps_interp(ji,boundary+1:nyfin-boundary) 
    497            wghts = wghts - (1. / ( (nb_connection_pts*irafx) - 1. ) )                       
    498         ENDDO 
    499         !                      
     499        DO ji = nxfin - boundary , nxfin - boundary - nb_connection_pts * irafx +1 , -1 
     500            wghts = wghts - (1. / (nb_connection_pts*irafx - 1. ) ) 
     501            DO jj=boundary+1,nyfin-boundary 
     502               wtab(ji,jj) = MAX(wghts, wtab(ji,jj)) 
     503            END DO 
     504        END DO   
     505 
    500506        wghts = 1. 
    501507        DO jj = boundary + 1 , boundary + nb_connection_pts * irafy 
    502            G1%gdepw_ps(boundary + nb_connection_pts * irafx + 1: & 
    503                 nxfin - boundary - nb_connection_pts * irafx ,jj) =          & 
    504                 (1. - wghts) * G1%gdepw_ps(boundary + nb_connection_pts * irafx + 1:  & 
    505                 nxfin - boundary - nb_connection_pts * irafx,jj) +                          & 
    506                 wghts * gdepw_ps_interp(boundary + nb_connection_pts * irafx + 1:   & 
    507                 nxfin - boundary - nb_connection_pts * irafx,jj) 
    508            wghts = wghts - (1. / (nb_connection_pts*irafx - 1. ) )                       
    509         ENDDO 
    510         !                     
     508            wghts = wghts - (1. / (nb_connection_pts*irafy - 1. ) ) 
     509            DO ji=boundary+1,nxfin-boundary 
     510               wtab(ji,jj) = MAX(wghts, wtab(ji,jj)) 
     511            END DO 
     512        END DO 
     513 
    511514        wghts = 1. 
    512         DO jj = nyfin - boundary , nyfin - boundary - nb_connection_pts * irafy+ 1 , -1 
    513            G1%gdepw_ps(boundary + nb_connection_pts * irafx + 1: & 
    514                 nxfin - boundary - nb_connection_pts * irafx ,jj) =                      & 
    515                 (1. - wghts) * G1%gdepw_ps(boundary + nb_connection_pts * irafx + 1: & 
    516                 nxfin - boundary - nb_connection_pts * irafx,jj) +                          & 
    517                 wghts * gdepw_ps_interp(boundary + nb_connection_pts * irafx + 1: & 
    518                 nxfin - boundary - nb_connection_pts * irafx,jj) 
    519            wghts = wghts - (1. / (nb_connection_pts*irafx - 1. ) )    
    520         ENDDO 
     515        DO jj = nyfin - boundary , nyfin - boundary - nb_connection_pts * irafy +1, -1 
     516            wghts = wghts - (1. / (nb_connection_pts*irafy - 1. ) ) 
     517            DO ji=boundary+1,nxfin-boundary 
     518               wtab(ji,jj) = MAX(wghts, wtab(ji,jj)) 
     519            END DO 
     520        END DO 
     521 
     522        G1%gdepw_ps(:,:) = (1.-wtab(:,:)) * G1%gdepw_ps(:,:) + gdepw_ps_interp(:,:)*wtab(:,:) 
    521523 
    522524        G1%bathy_meter = G1%gdepw_ps 
    523525        !                      
    524         connectionsize = 3 
     526        connectionsize = 2  
    525527        !  
    526528        IF(smoothing) THEN  
  • branches/2017/dev_merge_2017/NEMOGCM/TOOLS/NESTING/src/agrif_partial_steps.f90

    r5656 r9149  
    127127    ENDIF 
    128128    gdepw(1) = 0.0   
     129    IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]    
     130       !                      
     131       DO i = 1, N-1 
     132          e3t(i) = gdepw(i+1)-gdepw(i) 
     133       END DO 
     134       e3t(N) = e3t(N-1) 
     135     
     136       DO i = 2, N 
     137          e3w(i) = gdept(i) - gdept(i-1) 
     138       END DO 
     139       e3w(1  ) = 2. * (gdept(1) - gdepw(1)) 
     140    END IF 
    129141    ! 
    130142    ! Initialization of constant 
     
    329341    ENDIF 
    330342    gdepw(1) = 0.0 
     343    IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]    
     344       !                      
     345       DO i = 1, N-1 
     346          e3t(i) = gdepw(i+1)-gdepw(i) 
     347       END DO 
     348       e3t(N) = e3t(N-1) 
     349    END IF 
    331350    ! 
    332351    diff = 0       
     
    508527             x = ii-i - xdiff/2. 
    509528             val = parentgrid%gdepw_ps(ipt,jpt)+slopex * x 
     529!! chanut: uncomment this to get nearest neighbor interpolation 
     530!!             val = parentgrid%gdepw_ps(ipt,jpt)           
    510531             gdepwtemp(ii,j) = val 
    511532             IF (gdepwtemp(ii,j) < mindepth) THEN 
     
    565586             y = jj-j - xdiff/2. 
    566587             val = gdepwtemp(i,j) + slopey*y 
     588!! chanut: uncomment this to get nearest neighbor interpolation 
     589!!             val = gdepwtemp(i,j) 
    567590             gdepwtemp(i,jj) = val      
    568591          ENDDO 
     
    777800    ENDIF          
    778801    !          
     802    gdepw(1)=0. 
     803    IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]    
     804       !                      
     805       DO i = 1, jpk-1 
     806          e3t(i) = gdepw(i+1)-gdepw(i) 
     807       END DO 
     808       e3t(jpk) = e3t(jpk-1) 
     809    END IF 
    779810    !                 
    780811    DO i = 1,jpk 
  • branches/2017/dev_merge_2017/NEMOGCM/TOOLS/NESTING/src/agrif_readwrite.f90

    r8656 r9149  
    269269    LOGICAL,OPTIONAL :: Pacifique 
    270270    TYPE(Coordinates) :: CoarseGrid,ChildGrid 
     271    REAL*8 :: zdel 
     272    zdel = 0.5 ! Offset in degrees to extend extraction of bathymetry data 
    271273    ! 
    272274    IF( Dims_Existence('lon',name) .AND. Dims_Existence('lat',name) ) THEN 
     
    307309       END WHERE 
    308310       !           
    309        i_min = MAXLOC(topo_lon,mask = topo_lon < MINVAL(ChildGrid%nav_lon)) 
    310        i_max = MINLOC(topo_lon,mask = topo_lon > MAXVAL(ChildGrid%nav_lon))                     
    311        j_min = MAXLOC(topo_lat,mask = topo_lat < MINVAL(ChildGrid%nav_lat)) 
    312        j_max = MINLOC(topo_lat,mask = topo_lat > MAXVAL(ChildGrid%nav_lat)) 
     311       i_min = MAXLOC(topo_lon,mask = topo_lon < MINVAL(ChildGrid%nav_lon)-zdel) 
     312       i_max = MINLOC(topo_lon,mask = topo_lon > MAXVAL(ChildGrid%nav_lon)+zdel)                     
     313       j_min = MAXLOC(topo_lat,mask = topo_lat < MINVAL(ChildGrid%nav_lat)-zdel) 
     314       j_max = MINLOC(topo_lat,mask = topo_lat > MAXVAL(ChildGrid%nav_lat)+zdel) 
    313315       !           
    314316       tabdim1 = ( SIZE(topo_lon) - i_min(1) + 1 ) + i_max(1)                     
     
    336338       status = nf90_open(name,NF90_NOWRITE,ncid) 
    337339       status = nf90_inq_varid(ncid,elevation_name,varid) 
     340       ! 
     341       IF (status/=nf90_noerr) THEN 
     342          WRITE(*,*)"Can't find variable: ", elevation_name 
     343          STOP 
     344       ENDIF 
    338345       !           
    339346       status=nf90_get_var(ncid,varid,CoarseGrid%Bathy_meter(1:SIZE(topo_lon)-i_min(1)+1,:), & 
     
    345352    ELSE 
    346353       ! 
    347        i_min = MAXLOC(topo_lon,mask = topo_lon < MINVAL(ChildGrid%nav_lon)) 
    348        i_max = MINLOC(topo_lon,mask = topo_lon > MAXVAL(ChildGrid%nav_lon)) 
    349        j_min = MAXLOC(topo_lat,mask = topo_lat < MINVAL(ChildGrid%nav_lat)) 
    350        j_max = MINLOC(topo_lat,mask = topo_lat > MAXVAL(ChildGrid%nav_lat)) 
     354       i_min = MAXLOC(topo_lon,mask = topo_lon < MINVAL(ChildGrid%nav_lon)-zdel) 
     355       i_max = MINLOC(topo_lon,mask = topo_lon > MAXVAL(ChildGrid%nav_lon)+zdel) 
     356       j_min = MAXLOC(topo_lat,mask = topo_lat < MINVAL(ChildGrid%nav_lat)-zdel) 
     357       j_max = MINLOC(topo_lat,mask = topo_lat > MAXVAL(ChildGrid%nav_lat)+zdel) 
    351358       !       
    352359       IF(i_min(1)-2 >= 1 .AND. i_max(1)+3 <= SIZE(topo_lon,1) ) THEN 
     
    380387       status = nf90_open(name,NF90_NOWRITE,ncid) 
    381388       status = nf90_inq_varid(ncid,elevation_name,varid) 
     389       IF (status/=nf90_noerr) THEN 
     390          WRITE(*,*)"Can't find variable: ", elevation_name 
     391          STOP 
     392       ENDIF 
     393 
    382394       status = nf90_get_var(ncid,varid,CoarseGrid%Bathy_meter, & 
    383395          &                  start=(/i_min(1),j_min(1)/),count=(/tabdim1,tabdim2/)) 
  • branches/2017/dev_merge_2017/NEMOGCM/TOOLS/NESTING/src/agrif_types.f90

    r6819 r9149  
    5555  REAL*8 psur,pa0,pa1,pa2,adatrj 
    5656  !        
    57   LOGICAL ldbletanh 
     57  LOGICAL ldbletanh,ln_e3_dep 
    5858  LOGICAL partial_steps,smoothing,bathy_update 
    5959  LOGICAL new_topo,removeclosedseas,dimg,iom_activated 
     
    7676  NAMELIST /nesting/imin,imax,jmin,jmax,rho,rhot,bathy_update,updated_parent_file       
    7777  ! 
    78   NAMELIST /vertical_grid/ppkth,ppacr,ppdzmin,pphmax,psur,pa0,pa1,N,ldbletanh,pa2,ppkth2,ppacr2 
     78  NAMELIST /vertical_grid/ppkth,ppacr,ppdzmin,pphmax,psur,pa0,pa1,N,ldbletanh,ln_e3_dep,pa2,ppkth2,ppacr2 
    7979  !  
    8080  NAMELIST /partial_cells/partial_steps,parent_bathy_meter,parent_batmet_name,e3zps_min,e3zps_rat       
     
    9090  NAMELIST /restart_trc/ restart_trc_file,interp_type  
    9191 
    92   INTEGER :: connectionsize = 3 
     92  INTEGER :: connectionsize = 2  
    9393  ! 
    9494CONTAINS 
Note: See TracChangeset for help on using the changeset viewer.