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 14013 for NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_top_interp.F90 – NEMO

Ignore:
Timestamp:
2020-12-02T16:23:27+01:00 (3 years ago)
Author:
jchanut
Message:

Manually merge with Christian's branch agrif top sponge/interp #2129

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_top_interp.F90

    r13352 r14013  
    4343      Agrif_SpecialValue    = 0._wp 
    4444      Agrif_UseSpecialValue = .TRUE. 
    45       l_vremap = ln_vert_remap 
     45      l_vremap              = ln_vert_remap 
    4646      ! 
    4747      CALL Agrif_Bc_variable( trn_id, procname=interptrn ) 
     
    6464      INTEGER  :: item 
    6565      ! vertical interpolation: 
    66       REAL(wp) :: zhtot 
    67       REAL(wp), DIMENSION(k1:k2,1:jptra) :: tabin 
    68       REAL(wp), DIMENSION(k1:k2) :: h_in, z_in 
     66      REAL(wp) :: zhtot, zwgt 
     67      REAL(wp), DIMENSION(k1:k2,1:jptra) :: tabin, tabin_i 
     68      REAL(wp), DIMENSION(k1:k2) :: z_in, h_in_i, z_in_i 
    6969      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out 
    7070      !!---------------------------------------------------------------------- 
     
    8585         END DO 
    8686 
    87          IF( l_vremap .OR. l_ini_child) THEN 
    88             ! Interpolate thicknesses 
     87         IF( l_vremap .OR. l_ini_child .OR. ln_zps ) THEN 
     88            ! Fill cell depths (i.e. gdept) to be interpolated 
    8989            ! Warning: these are masked, hence extrapolated prior interpolation. 
    90             DO jk=k1,k2 
    91                DO jj=j1,j2 
    92                   DO ji=i1,i2 
    93                       ptab(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 
    94  
     90            DO jj=j1,j2 
     91               DO ji=i1,i2 
     92                  ptab(ji,jj,k1,jptra+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kmm_a) 
     93                  DO jk=k1+1,k2 
     94                     ptab(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * & 
     95                        & ( ptab(ji,jj,jk-1,jptra+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kmm_a)+e3t(ji,jj,jk,Kmm_a)) ) 
    9596                  END DO 
    9697               END DO 
    9798            END DO 
    9899 
    99             ! Extrapolate thicknesses in partial bottom cells: 
    100             ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
    101             IF (ln_zps) THEN 
    102                DO jj=j1,j2 
    103                   DO ji=i1,i2 
    104                       jk = mbkt(ji,jj) 
    105                       ptab(ji,jj,jk,jptra+1) = 0._wp 
    106                   END DO 
    107                END DO            
    108             END IF 
    109          
    110100            ! Save ssh at last level: 
    111101            IF (.NOT.ln_linssh) THEN 
    112102               ptab(i1:i2,j1:j2,k2,jptra+1) = ssh(i1:i2,j1:j2,Kmm_a)*tmask(i1:i2,j1:j2,1)  
    113             ELSE 
    114                ptab(i1:i2,j1:j2,k2,jptra+1) = 0._wp 
    115103            END IF       
    116104         ENDIF 
     
    126114            DO jj=j1,j2 
    127115               DO ji=i1,i2 
    128                   tr(ji,jj,:,:,Krhs_a) = 0.                   
     116                  tr(ji,jj,:,:,Krhs_a) = 0.   
     117                  ! 
     118                  ! Build vertical grids: 
    129119                  N_in = mbkt_parent(ji,jj) 
    130                   zhtot = 0._wp 
    131                   DO jk=1,N_in !k2 = jpk of parent grid 
    132                      IF (jk==N_in) THEN 
    133                         h_in(jk) = ht0_parent(ji,jj) + ptab(ji,jj,k2,n2) - zhtot 
    134                      ELSE 
    135                         h_in(jk) = ptab(ji,jj,jk,n2) 
    136                      ENDIF 
    137                      zhtot = zhtot + h_in(jk) 
    138                      tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1) 
    139                   END DO 
    140                   z_in(1) = 0.5_wp * h_in(1) - zhtot + ht0_parent(ji,jj) 
     120                  ! Input grid (account for partial cells if any): 
     121                  DO jk=1,N_in 
     122                     z_in(jk) = ptab(ji,jj,jk,n2) - ptab(ji,jj,k2,n2) 
     123                     tabin(jk,1:jptra) = ptab(ji,jj,jk,1:jptra) 
     124                  END DO 
     125                   
     126                  ! Intermediate grid: 
     127                  DO jk = 1, N_in 
     128                     h_in_i(jk) = e3t0_parent(ji,jj,jk) * &  
     129                       &       (1._wp + ptab(ji,jj,k2,n2)/(ht0_parent(ji,jj)*ssmask(ji,jj) + 1._wp - ssmask(ji,jj))) 
     130                  END DO 
     131                  z_in_i(1) = 0.5_wp * h_in_i(1) 
    141132                  DO jk=2,N_in 
    142                      z_in(jk) = z_in(jk-1) + 0.5_wp * (h_in(jk-1)+h_in(jk)) 
    143                   END DO 
    144  
    145                   N_out = 0 
    146                   DO jk=1,jpk ! jpk of child grid 
    147                      IF (tmask(ji,jj,jk) == 0._wp) EXIT  
    148                      N_out = N_out + 1 
     133                     z_in_i(jk) = z_in_i(jk-1) + 0.5_wp * ( h_in_i(jk) + h_in_i(jk-1) ) 
     134                  END DO 
     135                  z_in_i(1:N_in) = z_in_i(1:N_in)  - ptab(ji,jj,k2,n2) 
     136 
     137                  ! Output (Child) grid: 
     138                  N_out = mbkt(ji,jj) 
     139                  DO jk=1,N_out 
    149140                     h_out(jk) = e3t(ji,jj,jk,Krhs_a) 
    150141                  END DO 
    151  
    152                   z_out(1) = 0.5_wp * h_out(1) - SUM(h_out(1:N_out)) + ht_0(ji,jj) 
     142                  z_out(1) = 0.5_wp * h_out(1) 
    153143                  DO jk=2,N_out 
    154                      z_out(jk) = z_out(jk-1) + 0.5_wp * (h_out(jk-1)+h_out(jk)) 
    155                   END DO 
     144                     z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) ) 
     145                  END DO 
     146                  IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Krhs_a)                
    156147 
    157148                  IF (N_in*N_out > 0) THEN 
     
    159150                        CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),tr(ji,jj,1:N_out,1:jptra,Krhs_a),          & 
    160151                                      &   z_out(1:N_out),N_in,N_out,jptra)   
    161                      ELSE  
    162                         CALL reconstructandremap(tabin(1:N_in,1:jptra),h_in(1:N_in),tr(ji,jj,1:N_out,1:jptra,Krhs_a),   & 
    163                                       &   h_out(1:N_out),N_in,N_out,jptra)   
     152                     ELSE   
     153                        CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),tabin_i(1:N_in,1:jptra),                     & 
     154                                     &   z_in_i(1:N_in),N_in,N_in,jptra) 
     155                        CALL reconstructandremap(tabin_i(1:N_in,1:jptra),h_in_i(1:N_in),tr(ji,jj,1:N_out,1:jptra,Krhs_a), & 
     156                                      &   h_out(1:N_out),N_in,N_out,jptra)    
    164157                     ENDIF 
    165158                  ENDIF 
     
    169162  
    170163         ELSE 
    171           
     164 
     165            IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells  
     166                                             ! linear vertical interpolation 
     167               DO jj=j1,j2 
     168                  DO ji=i1,i2 
     169                     ! 
     170                     N_in  = mbkt(ji,jj) 
     171                     N_out = mbkt(ji,jj) 
     172                     z_in(1) = ptab(ji,jj,1,n2) 
     173                     tabin(1,1:jptra) = ptab(ji,jj,1,1:jptra) 
     174                     DO jk=2, N_in 
     175                        z_in(jk) = ptab(ji,jj,jk,n2) 
     176                        tabin(jk,1:jptra) = ptab(ji,jj,jk,1:jptra) 
     177                     END DO 
     178                     IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - ptab(ji,jj,k2,n2) 
     179                     z_out(1) = 0.5_wp * e3t(ji,jj,1,Krhs_a) 
     180                     DO jk=2, N_out 
     181                        z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Krhs_a) + e3t(ji,jj,jk,Krhs_a)) 
     182                     END DO 
     183                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Krhs_a) 
     184                     CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),ptab(ji,jj,1:N_out,1:jptra), & 
     185                                   &   z_out(1:N_out),N_in,N_out,jptra)   
     186                  END DO 
     187               END DO 
     188 
     189            ENDIF 
     190 
    172191            DO jn=1, jptra 
    173192                tr(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk)  
Note: See TracChangeset for help on using the changeset viewer.