Changeset 8861


Ignore:
Timestamp:
2017-11-30T16:58:38+01:00 (3 years ago)
Author:
timgraham
Message:

Added vertical interpolation of avm

Location:
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90

    r8822 r8861  
    576576      !!----------------------------------------------------------------------   
    577577      ! 
    578       return 
    579578       
    580579      zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp ) 
     
    602601      ! 
    603602      INTEGER  ::   ji, jj, jk, jn, iref, jref   ! dummy loop indices 
    604       INTEGER  ::   imin, imax, jmin, jmax 
     603      INTEGER  ::   imin, imax, jmin, jmax, N_in, N_out 
    605604      REAL(wp) ::   zrhox , zalpha1, zalpha2, zalpha3 
    606605      REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7 
    607606      LOGICAL :: western_side, eastern_side,northern_side,southern_side 
    608 ! VERTICAL REFINEMENT BEGIN 
    609607      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child 
    610608      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 
    611609      REAL(wp) :: h_in(k1:k2) 
    612610      REAL(wp) :: h_out(1:jpk) 
    613       INTEGER :: N_in, N_out 
    614       REAL(wp) :: h_diff 
    615       REAL(wp) :: zrhoxy 
    616 ! VERTICAL REFINEMENT END 
     611      REAL(wp) :: h_diff, zrhoxy 
    617612 
    618613      zrhoxy = Agrif_rhox()*Agrif_rhoy() 
    619  
    620  
    621  
    622614      IF (before) THEN          
    623             IF(Agrif_UseSpecialValue) THEN  
    624 !               Agrif_SpecialValue = -999._wp 
    625                Agrif_SpecialValue = 0._wp 
    626             ELSE 
    627                Agrif_SpecialValue = 0._wp 
    628             ENDIF 
    629             DO jn = n1,n2-1 
    630                DO jk=k1,k2 
    631                   DO jj=j1,j2 
    632                      DO ji=i1,i2 
    633                         ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk)! * tmask(ji,jj,jk) - & 
    634 !                                          & (tmask(ji,jj,jk)-1) * Agrif_SpecialValue 
    635                      END DO 
    636                   END DO 
    637                END DO 
    638             END DO 
     615         DO jn = n1,n2-1 
    639616            DO jk=k1,k2 
    640617               DO jj=j1,j2 
    641                   DO ji=i1,i2 
    642                      ptab(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk)  
    643                   END DO 
    644                END DO 
    645             END DO 
     618                 DO ji=i1,i2 
     619                       ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 
     620                 END DO 
     621              END DO 
     622           END DO 
     623        END DO 
     624        DO jk=k1,k2 
     625           DO jj=j1,j2 
     626              DO ji=i1,i2 
     627                 ptab(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk)  
     628              END DO 
     629           END DO 
     630        END DO 
    646631      ELSE  
    647632         western_side  = (nb == 1).AND.(ndir == 1) 
     
    649634         southern_side = (nb == 2).AND.(ndir == 1) 
    650635         northern_side = (nb == 2).AND.(ndir == 2) 
    651          Agrif_SpecialValue = 0._wp !reset now interpolation is done 
    652 ! VERTICAL REFINEMENT BEGIN 
    653636#ifdef key_vertical          
    654637         ptab_child(:,:,:,:) = 0. 
    655          do jj=j1,j2 
    656          do ji=i1,i2 
    657            iref = ji 
    658            jref = jj 
    659            if(western_side) iref=2 
    660            if(eastern_side) iref=nlci-1 
    661            if(southern_side) jref=2 
    662            if(northern_side) jref=nlcj-1 
    663            N_in = 0 
    664            DO jk=k1,k2 !k2 = jpk of parent grid 
    665              IF (ptab(ji,jj,jk,n2) == 0) EXIT 
    666              N_in = N_in + 1 
    667              tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)/(ptab(ji,jj,jk,n2)) 
    668              h_in(N_in) = ptab(ji,jj,jk,n2)/(e1e2t(ji,jj)*zrhoxy) 
    669            END DO 
    670            N_out = 0 
    671            DO jk=1,jpk ! jpk of child grid 
    672              IF (tmask(iref,jref,jk) == 0) EXIT ! TODO: Will not work with ISF. !This doesn't seem to work at the moment in GYRE but is OK in overflow model 
    673              N_out = N_out + 1 
    674              h_out(jk) = e3t_n(iref,jref,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 
    675            ENDDO 
    676            IF (N_in > 0) THEN 
    677              h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    678 !             if (h_diff > 0) then 
    679 !               h_in(N_in+1) = h_diff 
    680 !               N_in = N_in + 1 
    681 !             else 
    682 !               h_out(N_out+1) = -h_diff 
    683 !               N_out = N_out + 1 
    684 !             endif  
    685            ptab(ji,jj,k2,:) = ptab(ji,jj,k2-1,:) !what is this line for????? 
    686            do jn=1,jpts 
    687              call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 
    688            enddo 
    689 !           if (abs(h_diff) > 1000.) then 
    690 !           do jn=1,jpts 
    691 !             do jk=1,N_out 
    692 !             print *,'AVANT APRES = ',ji,jj,jk,N_out,ptab(ji,jj,jk,jn),ptab_child(ji,jj,jk,jn) 
    693 !             enddo 
    694 !           enddo 
    695 !         endif 
    696            ENDIF 
    697          enddo 
    698          enddo 
     638         DO jj=j1,j2 
     639            DO ji=i1,i2 
     640               iref = ji 
     641               jref = jj 
     642               if(western_side) iref=2 
     643               if(eastern_side) iref=nlci-1 
     644               if(southern_side) jref=2 
     645               if(northern_side) jref=nlcj-1 
     646               N_in = 0 
     647               DO jk=k1,k2 !k2 = jpk of parent grid 
     648                  IF (ptab(ji,jj,jk,n2) == 0) EXIT 
     649                  N_in = N_in + 1 
     650                  tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)/(ptab(ji,jj,jk,n2)) 
     651                  h_in(N_in) = ptab(ji,jj,jk,n2)/(e1e2t(ji,jj)*zrhoxy) 
     652               END DO 
     653               N_out = 0 
     654               DO jk=1,jpk ! jpk of child grid 
     655                  IF (tmask(iref,jref,jk) == 0) EXIT  
     656                  N_out = N_out + 1 
     657                  h_out(jk) = e3t_n(iref,jref,jk) 
     658               ENDDO 
     659               IF (N_in > 0) THEN 
     660                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     661                  DO jn=1,jpts 
     662                     call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 
     663                  ENDDO 
     664               ENDIF 
     665            ENDDO 
     666         ENDDO 
    699667#else 
    700          do jk=k1,k2 
    701             do jj=j1,j2 
    702                do ji=i1,i2 
    703                  ! This will be slow - Is there a way to speed it up and avoid divide by zero? 
    704                  IF (ptab(ji,jj,jk,n2) .NE. 0) THEN 
    705                     ptab_child(ji,jj,jk,n1:n2-1) = ptab(ji,jj,jk,n1:n2-1)/ptab(ji,jj,jk,n2) 
    706                  ELSE 
    707                     ptab_child(ji,jj,jk,n1:n2-1) = 0._wp 
    708                  ENDIF 
    709                enddo 
    710             enddo 
    711          enddo 
     668         DO jk=k1,k2 
     669            DO jj=j1,j2 
     670               DO ji=i1,i2 
     671                  ptab_child(ji,jj,jk,n1:n2-1) = ptab(ji,jj,jk,n1:n2-1)/(zrhoxy*e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
     672               ENDDO 
     673            ENDDO 
     674         ENDDO 
    712675#endif 
    713     
    714  
    715 ! VERTICAL REFINEMENT END 
    716  
    717676         ! 
    718677         zrhox = Agrif_Rhox() 
     
    996955      INTEGER :: ji,jj,jk 
    997956      REAL(wp) :: zrhox 
    998 ! VERTICAL REFINEMENT BEGIN 
    999957      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child 
    1000958      REAL(wp), DIMENSION(k1:k2) :: tabin 
     
    1006964      INTEGER :: jref 
    1007965 
    1008 ! VERTICAL REFINEMENT END 
    1009966      !!---------------------------------------------     
    1010967      !       
     
    1029986         ptab_child(:,:,:) = 0. 
    1030987#ifdef key_vertical 
    1031 ! VERTICAL REFINEMENT BEGIN 
    1032988         southern_side = (nb == 2).AND.(ndir == 1) 
    1033989         northern_side = (nb == 2).AND.(ndir == 2) 
     
    10611017                 CYCLE 
    10621018               ENDIF 
    1063  
    1064 !               IF (N_in * N_out > 0) THEN 
    1065 !                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    1066 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 
    1067 !                  if (h_diff < 0.) then 
    1068 !                     print *,'CHECK YOUR BATHY interpvn...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 
    1069 !                     stop 
    1070 !                  endif 
    1071 !               ENDIF 
    10721019               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
    1073           
    10741020            enddo 
    10751021         enddo 
    1076 ! in the following 
    1077 ! remove division of va by fs e3v, zrhox and e1v (already done) 
    10781022            DO jk=1,jpkm1 
    10791023               DO jj=j1,j2 
     
    10811025               END DO 
    10821026            END DO 
    1083 ! VERTICAL REFINEMENT END 
    10841027#else  
    10851028            DO jk=1,jpkm1 
     
    14481391# if defined key_zdftke 
    14491392 
    1450    SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, before ) 
     1393   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before ) 
    14511394      !!---------------------------------------------------------------------- 
    14521395      !!                  ***  ROUTINE interavm  *** 
    14531396      !!----------------------------------------------------------------------   
    1454       INTEGER                              , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2 
    1455       REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab 
     1397      INTEGER                              , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, m1, m2 
     1398      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) ::   ptab 
    14561399      LOGICAL                              , INTENT(in   ) ::   before 
     1400      REAL(wp), DIMENSION(k1:k2) :: tabin 
     1401      REAL(wp) :: h_in(k1:k2) 
     1402      REAL(wp) :: h_out(1:jpk) 
     1403      REAL(wp) :: zrhoxy 
     1404      INTEGER  :: N_in, N_out, ji, jj, jk 
    14571405      !!----------------------------------------------------------------------   
    14581406      !       
    1459       IF( before ) THEN 
    1460          ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2) 
    1461       ELSE 
    1462          avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
     1407      zrhoxy = Agrif_rhox()*Agrif_rhoy() 
     1408      IF (before) THEN          
     1409         DO jk=k1,k2 
     1410            DO jj=j1,j2 
     1411              DO ji=i1,i2 
     1412                    ptab(ji,jj,jk,1) = avm_k(ji,jj,jk) 
     1413              END DO 
     1414           END DO 
     1415        END DO 
     1416#ifdef key_vertical          
     1417        DO jk=k1,k2 
     1418           DO jj=j1,j2 
     1419              DO ji=i1,i2 
     1420                 ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e1e2t(ji,jj) * e3w_n(ji,jj,jk)  
     1421              END DO 
     1422           END DO 
     1423        END DO 
     1424#else 
     1425      ptab(i1:i2,j1:j2,k1:k2,2) = 0._wp 
     1426#endif 
     1427      ELSE  
     1428#ifdef key_vertical          
     1429         avm_k(i1:i2,j1:j2,1:jpk) = 0. 
     1430         DO jj=j1,j2 
     1431            DO ji=i1,i2 
     1432               N_in = 0 
     1433               DO jk=k1,k2 !k2 = jpk of parent grid 
     1434                  IF (ptab(ji,jj,jk,2) == 0) EXIT 
     1435                  N_in = N_in + 1 
     1436                  tabin(jk) = ptab(ji,jj,jk,1) 
     1437                  h_in(N_in) = ptab(ji,jj,jk,2)/(e1e2t(ji,jj)*zrhoxy) 
     1438               END DO 
     1439               N_out = 0 
     1440               DO jk=1,jpk ! jpk of child grid 
     1441                  IF (wmask(ji,jj,jk) == 0) EXIT  
     1442                  N_out = N_out + 1 
     1443                  h_out(jk) = e3t_n(ji,jj,jk) 
     1444               ENDDO 
     1445               IF (N_in > 0) THEN 
     1446                  CALL reconstructandremap(tabin(1:N_in),h_in,avm_k(ji,jj,1:N_out),h_out,N_in,N_out) 
     1447               ENDIF 
     1448            ENDDO 
     1449         ENDDO 
     1450#else 
     1451         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1) 
     1452#endif 
    14631453      ENDIF 
    14641454      ! 
  • branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_user.F90

    r8835 r8861  
    374374   CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/), en_id) 
    375375   CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/),avt_id) 
    376    CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/),avm_id) 
     376   CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/jpi,jpj,jpk,2/),avm_id) 
    377377# endif 
    378378 
     
    757757   ENDIF 
    758758 
     759   AGRIF_Parent(ln_vertical)=ln_vertical ! Will this work? 
     760 
    759761   CALL Agrif_Update_trc(0) 
    760762   ! 
     
    832834   INTEGER  ::   ios                 ! Local integer output status for namelist read 
    833835   INTEGER  ::   iminspon 
    834    NAMELIST/namagrif/ nn_cln_update, rn_sponge_tra, rn_sponge_dyn, ln_spc_dyn, ln_chk_bathy 
     836   NAMELIST/namagrif/ nn_cln_update, rn_sponge_tra, rn_sponge_dyn, ln_spc_dyn, ln_chk_bathy, ln_vertical 
    835837   !!-------------------------------------------------------------------------------------- 
    836838   ! 
     
    854856      WRITE(numout,*) '      use special values for dynamics   ln_spc_dyn    = ', ln_spc_dyn 
    855857      WRITE(numout,*) '      check bathymetry                  ln_chk_bathy  = ', ln_chk_bathy 
     858      WRITE(numout,*) '      use vertical interpolation        ln_vertical   = ', ln_vertical 
    856859      WRITE(numout,*)  
    857860   ENDIF 
Note: See TracChangeset for help on using the changeset viewer.