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 8135 for branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90 – NEMO

Ignore:
Timestamp:
2017-06-05T11:01:20+02:00 (7 years ago)
Author:
timgraham
Message:

Added changes to code in NST_SRC from dev_r5803_UKMO_AGRIF_Vert_interp

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90

    r7646 r8135  
    158158      INTEGER, INTENT(in) :: kt 
    159159      !        
     160      return 
     161       
    160162      IF( ( Agrif_NbStepint() .NE. 0 ) .AND. (kt /= 0) ) RETURN 
    161163#  if defined TWO_WAY 
    … …  
    185187      !! 
    186188      INTEGER :: ji,jj,jk,jn 
    187       !!--------------------------------------------- 
    188       ! 
    189       IF (before) THEN 
    190          DO jn = n1,n2 
     189! VERTICAL REFINEMENT BEGIN 
     190      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child 
     191      REAL(wp) :: h_in(k1:k2) 
     192      REAL(wp) :: h_out(1:jpk) 
     193      INTEGER  :: N_in, N_out 
     194      REAL(wp) :: h_diff 
     195      REAL(wp) :: zrho_xy 
     196      REAL(wp) :: tabin(k1:k2,n1:n2) 
     197! VERTICAL REFINEMENT END 
     198      !!--------------------------------------------- 
     199      ! 
     200      IF (before) THEN 
     201         zrho_xy = Agrif_rhox() * Agrif_rhoy()  
     202         DO jn = n1,n2-1 
    191203            DO jk=k1,k2 
    192204               DO jj=j1,j2 
    193205                  DO ji=i1,i2 
    194                      tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) 
     206                     tabres(ji,jj,jk,jn) = zrho_xy * tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 
    195207                  END DO 
    196208               END DO 
    197209            END DO 
    198210         END DO 
    199       ELSE 
     211         DO jk=k1,k2 
     212            DO jj=j1,j2 
     213               DO ji=i1,i2 
     214                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * zrho_xy * e1e2t(ji,jj) * e3t_n(ji,jj,jk)  
     215               END DO 
     216            END DO 
     217         END DO 
     218          
     219      ELSE 
     220! VERTICAL REFINEMENT BEGIN 
     221         tabres_child(:,:,:,:) = 0. 
     222 
     223         DO jj=j1,j2 
     224         DO ji=i1,i2 
     225           N_in = 0 
     226           DO jk=k1,k2 !k2 = jpk of child grid 
     227             IF (tabres(ji,jj,jk,n2) == 0) EXIT 
     228             N_in = N_in + 1 
     229             tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
     230             h_in(N_in) = tabres(ji,jj,jk,n2)/e1e2t(ji,jj) 
     231           ENDDO 
     232           N_out = 0 
     233           DO jk=1,jpk ! jpk of parent grid 
     234             IF (tmask(ji,jj,jk) == 0) EXIT ! TODO: Will not work with ISF 
     235             N_out = N_out + 1 
     236             h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above 
     237           ENDDO 
     238           IF (N_in > 0) THEN 
     239             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     240! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 
     241!             IF abs(h_diff > 1.e-8) THEN 
     242!               N_in = N_in + 1 
     243!               h_in(N_in) = h_diff 
     244!               tabin(N_in,:) = tabin(N_in-1,:) 
     245             IF (h_diff < 0) THEN 
     246             print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out)) 
     247             print *,'Nval = ',N_out,mbathy(ji,jj) 
     248             print *,'BATHY = ',gdepw_0(ji,jj,mbathy(ji,jj)+1),sum(e3t_0(ji,jj,1:mbathy(ji,jj))) 
     249             STOP 
     250!               N_out = N_out + 1 
     251!               h_out(N_out) = - h_diff 
     252             ENDIF 
     253             DO jn=n1,n2-1 
     254               CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 
     255             ENDDO 
     256          ENDIF 
     257         ENDDO 
     258         ENDDO 
     259          
     260! WARNING : 
     261! tabres replaced by tabres_child in the following 
     262! k1 k2 replaced by 1 jpk 
     263! VERTICAL REFINEMENT END 
     264 
    200265         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    201266            ! Add asselin part 
    202             DO jn = n1,n2 
    203                DO jk=k1,k2 
     267            DO jn = n1,n2-1 
     268               DO jk=1,jpk 
    204269                  DO jj=j1,j2 
    205270                     DO ji=i1,i2 
    206                         IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 
     271                        IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 
    207272                           tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) &  
    208                                  & + atfp * ( tabres(ji,jj,jk,jn) & 
     273                                 & + atfp * ( tabres_child(ji,jj,jk,jn) & 
    209274                                 &             - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
    210275                        ENDIF 
    … …  
    214279            ENDDO 
    215280         ENDIF 
    216          DO jn = n1,n2 
    217             DO jk=k1,k2 
     281         DO jn = n1,n2-1 
     282            DO jk=1,jpk 
    218283               DO jj=j1,j2 
    219284                  DO ji=i1,i2 
    220                      IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN  
    221                         tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) * tmask(ji,jj,jk) 
     285                     IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN  
     286                        tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk) 
    222287                     END IF 
    223288                  END DO 
    … …  
    230295 
    231296 
    232    SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, before ) 
     297   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    233298      !!--------------------------------------------- 
    234299      !!           *** ROUTINE updateu *** 
    235300      !!--------------------------------------------- 
    236       INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2 
    237       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres 
    238       LOGICAL                               , INTENT(in   ) :: before 
     301      INTEGER                                 , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
     302      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,2), INTENT(inout) :: tabres 
     303      LOGICAL                                 , INTENT(in   ) :: before 
    239304      ! 
    240305      INTEGER  ::   ji, jj, jk 
    241306      REAL(wp) ::   zrhoy 
     307! VERTICAL REFINEMENT BEGIN 
     308      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 
     309      REAL(wp) :: h_in(k1:k2) 
     310      REAL(wp) :: h_out(1:jpk) 
     311      INTEGER :: N_in, N_out 
     312      REAL(wp) :: h_diff 
     313      REAL(wp) :: tabin(k1:k2) 
     314! VERTICAL REFINEMENT END 
    242315      !!--------------------------------------------- 
    243316      !  
    244317      IF( before ) THEN 
    245318         zrhoy = Agrif_Rhoy() 
    246          DO jk = k1, k2 
    247             tabres(i1:i2,j1:j2,jk) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk) 
    248          END DO 
    249       ELSE 
    250319         DO jk=k1,k2 
    251320            DO jj=j1,j2 
    252321               DO ji=i1,i2 
    253                   tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk) 
     322                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk) 
     323                  tabres(ji,jj,jk,2) = umask(ji,jj,jk) * zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk)  
     324               END DO 
     325            END DO 
     326         END DO 
     327      ELSE 
     328! VERTICAL REFINEMENT BEGIN 
     329         tabres_child(:,:,:) = 0. 
     330          
     331         DO jj=j1,j2 
     332         DO ji=i1,i2 
     333           N_in = 0 
     334           DO jk=k1,k2 !k2=jpk of child grid 
     335             IF (tabres(ji,jj,jk,2) == 0) EXIT 
     336             N_in = N_in + 1 
     337             tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     338             h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj) 
     339           ENDDO 
     340           N_out = 0 
     341           DO jk=1,jpk 
     342             IF (umask(ji,jj,jk) == 0) EXIT 
     343             N_out = N_out + 1 
     344             h_out(N_out) = e3u_n(ji,jj,jk) 
     345           ENDDO 
     346           IF (N_in * N_out > 0) THEN 
     347             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     348! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 
     349             if (h_diff < 0.) then 
     350             print *,'CHECK YOUR BATHY ...' 
     351             stop 
     352!             else ! Extends with 0 
     353!             N_in = N_in + 1 
     354!             tabin(N_in) = 0. 
     355!             h_in(N_in) = h_diff 
     356             endif 
     357             CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
     358          ENDIF 
     359         ENDDO 
     360         ENDDO 
     361          
     362! WARNING : 
     363! tabres replaced by tabres_child in the following 
     364! k1 k2 replaced by 1 jpk 
     365! remove division by fs e3u (already done) 
     366! VERTICAL REFINEMENT END 
     367 
     368         DO jk=1,jpk 
     369            DO jj=j1,j2 
     370               DO ji=i1,i2 
     371!Following line now replaced by division higher up in vertical refinement case 
     372!                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk) 
    254373                  ! 
    255374                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    256375                     ub(ji,jj,jk) = ub(ji,jj,jk) &  
    257                            & + atfp * ( tabres(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 
     376                           & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 
    258377                  ENDIF 
    259378                  ! 
    260                   un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk) 
     379                  un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk) 
    261380               END DO 
    262381            END DO 
    … …  
    271390      !!           *** ROUTINE updatev *** 
    272391      !!--------------------------------------------- 
    273       INTEGER :: i1,i2,j1,j2,k1,k2 
     392      INTEGER :: i1,i2,j1,j2,k1,k2,n1,n2 
    274393      INTEGER :: ji,jj,jk 
    275       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: tabres 
     394      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,2) :: tabres 
    276395      LOGICAL :: before 
    277396      !! 
    278397      REAL(wp) :: zrhox 
     398! VERTICAL REFINEMENT BEGIN 
     399      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 
     400      REAL(wp) :: h_in(k1:k2) 
     401      REAL(wp) :: h_out(1:jpk) 
     402      INTEGER :: N_in, N_out 
     403      REAL(wp) :: h_diff 
     404      REAL(wp) :: tabin(k1:k2) 
     405! VERTICAL REFINEMENT END 
    279406      !!---------------------------------------------       
    280407      ! 
    … …  
    284411            DO jj=j1,j2 
    285412               DO ji=i1,i2 
    286                   tabres(ji,jj,jk) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) 
    287                END DO 
    288             END DO 
    289          END DO 
    290       ELSE 
    291          DO jk=k1,k2 
     413                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) 
     414                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk)  
     415               END DO 
     416            END DO 
     417         END DO 
     418      ELSE 
     419! VERTICAL REFINEMENT BEGIN 
     420         tabres_child(:,:,:) = 0. 
     421          
     422         DO jj=j1,j2 
     423         DO ji=i1,i2 
     424           N_in = 0 
     425           DO jk=k1,k2 
     426             IF (tabres(ji,jj,jk,2) == 0) EXIT 
     427             N_in = N_in + 1 
     428             tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     429             h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj) 
     430           ENDDO 
     431           N_out = 0 
     432           DO jk=1,jpk 
     433             IF (vmask(ji,jj,jk) == 0) EXIT 
     434             N_out = N_out + 1 
     435             h_out(N_out) = e3v_n(ji,jj,jk) 
     436           ENDDO 
     437           IF (N_in * N_out > 0) THEN 
     438             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     439! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 
     440             if (h_diff < 0.) then 
     441             print *,'CHECK YOUR BATHY ...' 
     442             stop 
     443!             else ! Extends with 0 
     444!             N_in = N_in + 1 
     445!             tabin(N_in) = 0. 
     446!             h_in(N_in) = h_diff 
     447             endif 
     448             CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
     449          ENDIF 
     450         ENDDO 
     451         ENDDO 
     452          
     453! WARNING : 
     454! tabres replaced by tabres_child in the following 
     455! k1 k2 replaced by 1 jpk 
     456! remove division by fs e3v (already done) 
     457! VERTICAL REFINEMENT END 
     458 
     459         DO jk=k1,jpk 
    292460            DO jj=j1,j2 
    293461               DO ji=i1,i2 
    294                   tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk) 
     462!Following line now replaced by division higher up in vertical refinement case 
     463!                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk) 
    295464                  ! 
    296465                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    297466                     vb(ji,jj,jk) = vb(ji,jj,jk) &  
    298                            & + atfp * ( tabres(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 
     467                           & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 
    299468                  ENDIF 
    300469                  ! 
    301                   vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk) 
     470                  vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk) 
    302471               END DO 
    303472            END DO 
Note: See TracChangeset for help on using the changeset viewer.