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 6258 for branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO – NEMO

Ignore:
Timestamp:
2016-01-15T13:11:56+01:00 (8 years ago)
Author:
timgraham
Message:

First inclusion of Laurent Debreu's modified code for vertical refinement.
Still a lot of outstanding issues:
1) conv preprocessor fails for limrhg.F90 at the moment (for now I've run without ice model)
2) conv preprocessor fails for STO code - removed this code from testing for now
3) conv preprocessor fails for cpl_oasis.F90 - can work round this by modifying code but the preprocessor should be fixed to deal with this.

After that code compiles and can be run for horizontal grid refinement. Not yet working for vertical refinement.

Location:
branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90

    r5819 r6258  
    1717 
    1818   PUBLIC agrif_oce_alloc ! routine called by nemo_init in nemogcm.F90 
     19   PUBLIC reconstructandremap 
    1920 
    2021   !                                              !!* Namelist namagrif: AGRIF parameters 
    … …  
    5859   INTEGER :: un_update_id, vn_update_id                              ! AGRIF profiles for udpates 
    5960   INTEGER :: tsn_sponge_id, un_sponge_id, vn_sponge_id               ! AGRIF profiles for sponge layers 
     61! VERTICAL REFINEMENT BEGIN 
     62   INTEGER :: scales_t_id, scales_u_id, scales_v_id 
     63! VERTICAL REFINEMENT END 
     64 
    6065# if defined key_top 
    6166   INTEGER :: trn_id, trn_sponge_id 
    … …  
    6469   INTEGER :: ub2b_update_id, vb2b_update_id 
    6570   INTEGER :: e3t_id, e1u_id, e2v_id, sshn_id 
    66    INTEGER :: scales_t_id 
    6771# if defined key_zdftke 
    6872   INTEGER :: avt_id, avm_id, en_id 
    … …  
    7377   !!---------------------------------------------------------------------- 
    7478   !! NEMO/NST 3.3.1 , NEMO Consortium (2011) 
    75    !! $Id$ 
     79   !! $Id: agrif_oce.F90 5081 2015-02-13 09:51:27Z smasson $ 
    7680   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    7781   !!---------------------------------------------------------------------- 
    … …  
    104108   END FUNCTION agrif_oce_alloc 
    105109 
     110      subroutine reconstructandremap(tabin,hin,tabout,hout,N,Nout)       
     111      implicit none 
     112      integer N, Nout 
     113      real tabin(N), tabout(Nout) 
     114      real hin(N), hout(Nout) 
     115      real coeffremap(N,3),zwork(N,3) 
     116      real zwork2(N+1,3) 
     117      integer k 
     118      double precision, parameter :: dsmll=1.0d-8   
     119      real q,q01,q02,q001,q002,q0 
     120      real z_win(1:N+1), z_wout(1:Nout+1) 
     121      real,parameter :: dpthin = 1.D-3 
     122      integer :: k1, kbox, ktop, ka, kbot 
     123      real :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop 
     124 
     125      z_win(1)=0.; z_wout(1)= 0. 
     126      do k=1,N 
     127       z_win(k+1)=z_win(k)+hin(k) 
     128      enddo  
     129       
     130      do k=1,Nout 
     131       z_wout(k+1)=z_wout(k)+hout(k)        
     132      enddo        
     133 
     134        do k=2,N 
     135          zwork(k,1)=1./(hin(k-1)+hin(k)) 
     136        enddo 
     137         
     138        do k=2,N-1 
     139          q0 = 1./(hin(k-1)+hin(k)+hin(k+1)) 
     140          zwork(k,2)=hin(k-1)+2.*hin(k)+hin(k+1) 
     141          zwork(k,3)=q0 
     142        enddo        
     143       
     144        do k= 2,N 
     145        zwork2(k,1)=zwork(k,1)*(tabin(k)-tabin(k-1)) 
     146        enddo 
     147         
     148        coeffremap(:,1) = tabin(:) 
     149  
     150         do k=2,N-1 
     151        q001 = hin(k)*zwork2(k+1,1) 
     152        q002 = hin(k)*zwork2(k,1)         
     153        if (q001*q002 < 0) then 
     154          q001 = 0. 
     155          q002 = 0. 
     156        endif 
     157        q=zwork(k,2) 
     158        q01=q*zwork2(k+1,1) 
     159        q02=q*zwork2(k,1) 
     160        if (abs(q001) > abs(q02)) q001 = q02 
     161        if (abs(q002) > abs(q01)) q002 = q01 
     162         
     163        q=(q001-q002)*zwork(k,3) 
     164        q001=q001-q*hin(k+1) 
     165        q002=q002+q*hin(k-1) 
     166         
     167        coeffremap(k,3)=coeffremap(k,1)+q001 
     168        coeffremap(k,2)=coeffremap(k,1)-q002 
     169         
     170        zwork2(k,1)=(2.*q001-q002)**2 
     171        zwork2(k,2)=(2.*q002-q001)**2 
     172        enddo 
     173         
     174        do k=1,N 
     175        if     (k.eq.1 .or. k.eq.N .or.   hin(k).le.dpthin) then 
     176        coeffremap(k,3) = coeffremap(k,1) 
     177        coeffremap(k,2) = coeffremap(k,1) 
     178        zwork2(k,1) = 0. 
     179        zwork2(k,2) = 0. 
     180        endif 
     181        enddo 
     182         
     183        do k=2,N 
     184        q002=max(zwork2(k-1,2),dsmll) 
     185        q001=max(zwork2(k,1),dsmll) 
     186        zwork2(k,3)=(q001*coeffremap(k-1,3)+q002*coeffremap(k,2))/(q001+q002) 
     187        enddo 
     188         
     189        zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3) 
     190        zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3) 
     191  
     192        do k=1,N 
     193        q01=zwork2(k+1,3)-coeffremap(k,1) 
     194        q02=coeffremap(k,1)-zwork2(k,3) 
     195        q001=2.*q01 
     196        q002=2.*q02 
     197        if (q01*q02<0) then 
     198          q01=0. 
     199          q02=0. 
     200        elseif (abs(q01)>abs(q002)) then 
     201          q01=q002 
     202        elseif (abs(q02)>abs(q001)) then 
     203          q02=q001 
     204        endif 
     205        coeffremap(k,2)=coeffremap(k,1)-q02 
     206        coeffremap(k,3)=coeffremap(k,1)+q01 
     207        enddo 
     208 
     209      zbot=0.0 
     210      kbot=1 
     211      do k=1,Nout 
     212        ztop=zbot  !top is bottom of previous layer 
     213        ktop=kbot 
     214        if     (ztop.ge.z_win(ktop+1)) then 
     215          ktop=ktop+1 
     216        endif 
     217         
     218        zbot=z_wout(k+1) 
     219        zthk=zbot-ztop 
     220 
     221        if     (zthk.gt.dpthin .and. ztop.lt.z_wout(Nout+1)) then 
     222 
     223          kbot=ktop 
     224          do while (z_win(kbot+1).lt.zbot.and.kbot.lt.N) 
     225            kbot=kbot+1 
     226          enddo 
     227          zbox=zbot 
     228          do k1= k+1,Nout 
     229            if     (z_wout(k1+1)-z_wout(k1).gt.dpthin) then 
     230              exit !thick layer 
     231            else 
     232              zbox=z_wout(k1+1)  !include thin adjacent layers 
     233              if     (zbox.eq.z_wout(Nout+1)) then 
     234                exit !at bottom 
     235              endif 
     236            endif 
     237          enddo 
     238          zthk=zbox-ztop 
     239 
     240          kbox=ktop 
     241          do while (z_win(kbox+1).lt.zbox.and.kbox.lt.N) 
     242            kbox=kbox+1 
     243          enddo 
     244           
     245          if     (ktop.eq.kbox) then 
     246 
     247 
     248            if     (z_wout(k)  .ne.z_win(kbox)   .or.z_wout(k+1).ne.z_win(kbox+1)     ) then 
     249 
     250              if     (hin(kbox).gt.dpthin) then 
     251                q001 = (zbox-z_win(kbox))/hin(kbox) 
     252                q002 = (ztop-z_win(kbox))/hin(kbox) 
     253                q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002) 
     254                q02=q01-1.+(q001+q002) 
     255                q0=1.-q01-q02 
     256              else 
     257                q0 = 1.0 
     258                q01 = 0. 
     259                q02 = 0. 
     260              endif 
     261          tabout(k)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3) 
     262               
     263            else 
     264            tabout(k) = tabin(kbox) 
     265               
     266            endif  
     267 
     268          else 
     269 
     270            if     (ktop.le.k .and. kbox.ge.k) then 
     271              ka = k 
     272            elseif (kbox-ktop.ge.3) then 
     273              ka = (kbox+ktop)/2 
     274            elseif (hin(ktop).ge.hin(kbox)) then 
     275              ka = ktop 
     276            else 
     277              ka = kbox 
     278            endif !choose ka 
     279 
     280            offset=coeffremap(ka,1) 
     281 
     282            qtop = z_win(ktop+1)-ztop !partial layer thickness 
     283            if     (hin(ktop).gt.dpthin) then 
     284              q=(ztop-z_win(ktop))/hin(ktop) 
     285              q01=q*(q-1.) 
     286              q02=q01+q 
     287              q0=1-q01-q02             
     288            else 
     289              q0 = 1. 
     290              q01 = 0. 
     291              q02 = 0. 
     292            endif 
     293             
     294            tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop 
     295 
     296            do k1= ktop+1,kbox-1 
     297              tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1) 
     298            enddo !k1 
     299 
     300             
     301            qbot = zbox-z_win(kbox) !partial layer thickness 
     302            if     (hin(kbox).gt.dpthin) then 
     303              q=qbot/hin(kbox) 
     304              q01=(q-1.)**2 
     305              q02=q01-1.+q 
     306              q0=1-q01-q02                             
     307            else 
     308              q0 = 1.0 
     309              q01 = 0. 
     310              q02 = 0. 
     311            endif 
     312            
     313            tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot 
     314             
     315            rpsum=1.0d0/zthk 
     316              tabout(k)=offset+tsum*rpsum 
     317               
     318          endif !single or multiple layers 
     319        else 
     320        if (k==1) then 
     321        write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(k+1),hout(1) 
     322        endif 
     323         tabout(k) = tabout(k-1) 
     324           
     325        endif !normal:thin layer 
     326      enddo !k 
     327             
     328      return 
     329      end subroutine reconstructandremap 
     330       
    106331#endif 
    107332   !!====================================================================== 

  • branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90

    r5819 r6258  
    3838 
    3939   PUBLIC   Agrif_tra, Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_ssh_ts, Agrif_dta_ts 
     40! VERTICAL REFINEMENT BEGIN    
     41   PUBLIC   Agrif_Init_InterpScales 
     42! VERTICAL REFINEMENT END 
    4043   PUBLIC   interpun, interpvn, interpun2d, interpvn2d  
    4144   PUBLIC   interptsn,  interpsshn 
    … …  
    5053   !!---------------------------------------------------------------------- 
    5154   !! NEMO/NST 3.6 , NEMO Consortium (2010) 
    52    !! $Id$ 
     55   !! $Id: agrif_opa_interp.F90 5081 2015-02-13 09:51:27Z smasson $ 
    5356   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    5457   !!---------------------------------------------------------------------- 
    5558 
     59! VERTICAL REFINEMENT BEGIN 
     60   REAL, DIMENSION(:,:,:), ALLOCATABLE :: interp_scales_t, interp_scales_u, interp_scales_v 
     61!$AGRIF_DO_NOT_TREAT 
     62   LOGICAL :: scaleT, scaleU, scaleV = .FALSE. 
     63!$AGRIF_END_DO_NOT_TREAT 
     64! VERTICAL REFINEMENT END 
     65 
    5666CONTAINS 
     67 
     68! VERTICAL REFINEMENT BEGIN 
     69 
     70   SUBROUTINE Agrif_Init_InterpScales 
     71 
     72    scaleT = .TRUE. 
     73    Call Agrif_Bc_Variable(scales_t_id,calledweight=1.,procname=interpscales) 
     74    scaleT = .FALSE. 
     75     
     76    scaleU = .TRUE. 
     77    Call Agrif_Bc_Variable(scales_u_id,calledweight=1.,procname=interpscales) 
     78    scaleU = .FALSE. 
     79 
     80    scaleV = .TRUE. 
     81    Call Agrif_Bc_Variable(scales_v_id,calledweight=1.,procname=interpscales) 
     82    scaleV = .FALSE. 
     83 
     84   END SUBROUTINE Agrif_Init_InterpScales 
     85    
     86   SUBROUTINE interpscales(ptab,i1,i2,j1,j2,k1,k2,before) 
     87      !!--------------------------------------------- 
     88      !!   *** ROUTINE interpscales *** 
     89      !!--------------------------------------------- 
     90       
     91      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 
     92      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 
     93 
     94      INTEGER :: ji, jj, jk 
     95      LOGICAL :: before 
     96 
     97      IF (before) THEN 
     98      IF (scaleT ) THEN 
     99      DO jk=k1,k2 
     100         DO jj=j1,j2 
     101            DO ji=i1,i2 
     102!               ptab(ji,jj,jk) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) 
     103               ptab(ji,jj,jk) = fse3t_n(ji,jj,jk) 
     104            END DO 
     105         END DO 
     106      END DO 
     107      ELSEIF (scaleU) THEN 
     108      DO jk=k1,k2 
     109         DO jj=j1,j2 
     110            DO ji=i1,i2 
     111!               ptab(ji,jj,jk) = fse3u_n(ji,jj,jk) * umask(ji,jj,jk) 
     112!               ptab(ji,jj,jk) = fse3u_n(ji,jj,jk) 
     113                ptab(ji,jj,jk) = umask(ji,jj,jk) 
     114            END DO 
     115         END DO 
     116      END DO 
     117      ELSEIF (scaleV) THEN 
     118      DO jk=k1,k2 
     119         DO jj=j1,j2 
     120            DO ji=i1,i2 
     121!               ptab(ji,jj,jk) = fse3v_n(ji,jj,jk) * vmask(ji,jj,jk) 
     122!               ptab(ji,jj,jk) = fse3v_n(ji,jj,jk) 
     123               ptab(ji,jj,jk) = vmask(ji,jj,jk) 
     124            END DO 
     125         END DO 
     126      END DO 
     127      ENDIF 
     128      ELSE 
     129      IF (scaleT ) THEN 
     130      IF (.not.allocated(interp_scales_t)) allocate(interp_scales_t(jpi,jpj,k1:k2)) 
     131      DO jk=k1,k2 
     132         DO jj=j1,j2 
     133            DO ji=i1,i2 
     134               interp_scales_t(ji,jj,jk) = ptab(ji,jj,jk) 
     135            END DO 
     136         END DO 
     137      END DO 
     138      ELSEIF (scaleU) THEN 
     139      IF (.not.allocated(interp_scales_u)) allocate(interp_scales_u(jpi,jpj,k1:k2)) 
     140      DO jk=k1,k2 
     141         DO jj=j1,j2 
     142            DO ji=i1,i2 
     143               interp_scales_u(ji,jj,jk) = ptab(ji,jj,jk) 
     144            END DO 
     145         END DO 
     146      END DO 
     147      ELSEIF (scaleV) THEN 
     148      IF (.not.allocated(interp_scales_v)) allocate(interp_scales_v(jpi,jpj,k1:k2)) 
     149      DO jk=k1,k2 
     150         DO jj=j1,j2 
     151            DO ji=i1,i2 
     152               interp_scales_v(ji,jj,jk) = ptab(ji,jj,jk) 
     153            END DO 
     154         END DO 
     155      END DO 
     156      ENDIF 
     157      ENDIF 
     158 
     159   END SUBROUTINE interpscales 
     160 
     161! VERTICAL REFINEMENT END 
    57162 
    58163   SUBROUTINE Agrif_tra 
    … …  
    611716      REAL(wp) ::   zalpha 
    612717      ! 
     718      return 
     719       
    613720      zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp ) 
    614721      IF( zalpha > 1. )   zalpha = 1. 
    … …  
    638745      REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7 
    639746      LOGICAL :: western_side, eastern_side,northern_side,southern_side 
     747! VERTICAL REFINEMENT BEGIN 
     748      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child 
     749      REAL(wp) :: h_in(k1:k2) 
     750      REAL(wp) :: h_out(1:jpk) 
     751      INTEGER :: N_in, N_out 
     752      REAL(wp) :: h_diff 
     753! VERTICAL REFINEMENT END 
    640754 
    641755      IF (before) THEN          
    642756         ptab(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2) 
    643       ELSE 
     757      ELSE  
     758! VERTICAL REFINEMENT BEGIN 
     759 
     760         ptab_child(:,:,:,:) = 0. 
     761         do jj=j1,j2 
     762         do ji=i1,i2 
     763           h_in(k1:k2) = interp_scales_t(ji,jj,k1:k2) 
     764           h_out(1:jpk) = fse3t(ji,jj,1:jpk) 
     765           h_diff = sum(h_out(1:jpk-1))-sum(h_in(k1:k2-1)) 
     766           N_in = k2-1 
     767           N_out = jpk-1 
     768           if (h_diff > 0) then 
     769             h_in(N_in+1) = h_diff 
     770             N_in = N_in + 1 
     771           else 
     772             h_out(N_out+1) = -h_diff 
     773             N_out = N_out + 1 
     774           endif  
     775           ptab(ji,jj,k2,:) = ptab(ji,jj,k2-1,:) 
     776           do jn=1,jpts 
     777             call reconstructandremap(ptab(ji,jj,1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 
     778           enddo 
     779!           if (abs(h_diff) > 1000.) then 
     780!           do jn=1,jpts 
     781!             do jk=1,N_out 
     782!             print *,'AVANT APRES = ',ji,jj,jk,N_out,ptab(ji,jj,jk,jn),ptab_child(ji,jj,jk,jn) 
     783!             enddo 
     784!           enddo 
     785!         endif 
     786         enddo 
     787         enddo 
     788 
     789! VERTICAL REFINEMENT END 
     790 
    644791         ! 
    645792         western_side  = (nb == 1).AND.(ndir == 1) 
    … …  
    671818         IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci-2         
    672819         ! 
     820! VERTICAL REFINEMENT BEGIN 
     821 
     822! WARNING : 
     823! ptab replaced by ptab_child in the following 
     824! k1 k2 replaced by 1 jpk 
     825! VERTICAL REFINEMENT END 
    673826         IF( eastern_side) THEN 
    674827            DO jn = 1, jpts 
    675                tsa(nlci,j1:j2,k1:k2,jn) = zalpha1 * ptab(nlci,j1:j2,k1:k2,jn) + zalpha2 * ptab(nlci-1,j1:j2,k1:k2,jn) 
     828               tsa(nlci,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(nlci,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(nlci-1,j1:j2,1:jpk,jn) 
    676829               DO jk = 1, jpkm1 
    677830                  DO jj = jmin,jmax 
    … …  
    692845         IF( northern_side ) THEN             
    693846            DO jn = 1, jpts 
    694                tsa(i1:i2,nlcj,k1:k2,jn) = zalpha1 * ptab(i1:i2,nlcj,k1:k2,jn) + zalpha2 * ptab(i1:i2,nlcj-1,k1:k2,jn) 
     847               tsa(i1:i2,nlcj,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,nlcj,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,nlcj-1,1:jpk,jn) 
    695848               DO jk = 1, jpkm1 
    696849                  DO ji = imin,imax 
    … …  
    711864         IF( western_side) THEN             
    712865            DO jn = 1, jpts 
    713                tsa(1,j1:j2,k1:k2,jn) = zalpha1 * ptab(1,j1:j2,k1:k2,jn) + zalpha2 * ptab(2,j1:j2,k1:k2,jn) 
     866               tsa(1,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(1,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(2,j1:j2,1:jpk,jn) 
    714867               DO jk = 1, jpkm1 
    715868                  DO jj = jmin,jmax 
    … …  
    729882         IF( southern_side ) THEN            
    730883            DO jn = 1, jpts 
    731                tsa(i1:i2,1,k1:k2,jn) = zalpha1 * ptab(i1:i2,1,k1:k2,jn) + zalpha2 * ptab(i1:i2,2,k1:k2,jn) 
     884               tsa(i1:i2,1,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,1,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,2,1:jpk,jn) 
    732885               DO jk=1,jpk       
    733886                  DO ji=imin,imax 
    … …  
    749902         ! East south 
    750903         IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 
    751             tsa(nlci-1,2,:,:) = ptab(nlci-1,2,:,:) 
     904            tsa(nlci-1,2,:,:) = ptab_child(nlci-1,2,:,:) 
    752905         ENDIF 
    753906         ! East north 
    754907         IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 
    755             tsa(nlci-1,nlcj-1,:,:) = ptab(nlci-1,nlcj-1,:,:) 
     908            tsa(nlci-1,nlcj-1,:,:) = ptab_child(nlci-1,nlcj-1,:,:) 
    756909         ENDIF 
    757910         ! West south 
    758911         IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 
    759             tsa(2,2,:,:) = ptab(2,2,:,:) 
     912            tsa(2,2,:,:) = ptab_child(2,2,:,:) 
    760913         ENDIF 
    761914         ! West north 
    762915         IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 
    763             tsa(2,nlcj-1,:,:) = ptab(2,nlcj-1,:,:) 
     916            tsa(2,nlcj-1,:,:) = ptab_child(2,nlcj-1,:,:) 
    764917         ENDIF 
    765918         ! 
    … …  
    794947   END SUBROUTINE interpsshn 
    795948 
    796    SUBROUTINE interpun(ptab,i1,i2,j1,j2,k1,k2, before) 
     949   SUBROUTINE interpun(ptab,i1,i2,j1,j2,k1,k2,m1,m2,before,nb,ndir) 
    797950      !!--------------------------------------------- 
    798951      !!   *** ROUTINE interpun *** 
    799952      !!---------------------------------------------     
    800953      !! 
    801       INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 
    802       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 
     954      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 
     955      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 
    803956      LOGICAL, INTENT(in) :: before 
     957      INTEGER, INTENT(in) :: nb , ndir 
    804958      !! 
    805959      INTEGER :: ji,jj,jk 
    806       REAL(wp) :: zrhoy  
     960      REAL(wp) :: zrhoy 
     961! VERTICAL REFINEMENT BEGIN 
     962      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child 
     963      REAL(wp), DIMENSION(k1:k2) :: tabin 
     964      REAL(wp) :: h_in(k1:k2) 
     965      REAL(wp) :: h_out(1:jpk) 
     966      INTEGER :: N_in, N_out 
     967      REAL(wp) :: h_diff 
     968      LOGICAL :: western_side, eastern_side 
     969      INTEGER :: iref 
     970 
     971! VERTICAL REFINEMENT END 
    807972      !!---------------------------------------------     
    808973      ! 
    … …  
    811976            DO jj=j1,j2 
    812977               DO ji=i1,i2 
    813                   ptab(ji,jj,jk) = e2u(ji,jj) * un(ji,jj,jk) 
    814                   ptab(ji,jj,jk) = ptab(ji,jj,jk) * fse3u(ji,jj,jk) 
     978                  ptab(ji,jj,jk,1) = e2u(ji,jj) * un(ji,jj,jk) 
     979                  ptab(ji,jj,jk,1) = ptab(ji,jj,jk,1) * fse3u(ji,jj,jk) 
     980                  ptab(ji,jj,jk,2) = fse3u(ji,jj,jk) 
    815981               END DO 
    816982            END DO 
    817983         END DO 
    818984      ELSE 
     985! VERTICAL REFINEMENT BEGIN 
     986         western_side  = (nb == 1).AND.(ndir == 1) 
     987         eastern_side  = (nb == 1).AND.(ndir == 2) 
     988          
     989         ptab_child(:,:,:) = 0. 
     990         do jj=j1,j2 
     991         do ji=i1,i2 
     992         iref = ji 
     993         IF (western_side) iref = 2 
     994         IF (eastern_side) iref = nlci-2 
     995 
     996         N_in = 0 
     997         do jk=k1,k2 
     998           if (interp_scales_u(ji,jj,jk) == 0) EXIT 
     999             N_in = N_in + 1 
     1000             tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 
     1001             h_in(N_in) = ptab(ji,jj,jk,2) 
     1002         enddo 
     1003          
     1004         IF (N_in == 0) THEN 
     1005           ptab_child(ji,jj,:) = 0. 
     1006           CYCLE 
     1007         ENDIF 
     1008          
     1009         N_out = 0 
     1010         do jk=1,jpk 
     1011           if (umask(iref,jj,jk) == 0) EXIT 
     1012           N_out = N_out + 1 
     1013           h_out(N_out) = fse3u(ji,jj,jk) 
     1014         enddo 
     1015          
     1016         IF (N_out == 0) THEN 
     1017           ptab_child(ji,jj,:) = 0. 
     1018           CYCLE 
     1019         ENDIF 
     1020          
     1021         h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     1022         IF (h_diff > 0.) THEN 
     1023           N_in = N_in + 1 
     1024           h_in(N_in) = h_diff 
     1025           tabin(N_in) = 0. 
     1026         ELSE 
     1027           h_out(N_out) = h_out(N_out) - h_diff 
     1028         ENDIF 
     1029          
     1030         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) 
     1031          
     1032         ptab_child(ji,jj,N_out) = ptab_child(ji,jj,N_out) * h_out(N_out) / fse3u(ji,jj,N_out) 
     1033 
     1034         ENDDO 
     1035         ENDDO 
     1036 
     1037! in the following 
     1038! remove division of ua by fs e3u (already done) 
     1039! VERTICAL REFINEMENT END 
     1040 
    8191041         zrhoy = Agrif_Rhoy() 
    8201042         DO jk=1,jpkm1 
    8211043            DO jj=j1,j2 
    822                ua(i1:i2,jj,jk) = (ptab(i1:i2,jj,jk)/(zrhoy*e2u(i1:i2,jj))) 
    823                ua(i1:i2,jj,jk) = ua(i1:i2,jj,jk) / fse3u(i1:i2,jj,jk) 
     1044               ua(i1:i2,jj,jk) = (ptab_child(i1:i2,jj,jk)/(zrhoy*e2u(i1:i2,jj))) 
    8241045            END DO 
    8251046         END DO 
    … …  
    8611082 
    8621083 
    863    SUBROUTINE interpvn(ptab,i1,i2,j1,j2,k1,k2, before) 
     1084   SUBROUTINE interpvn(ptab,i1,i2,j1,j2,k1,k2,m1,m2,before,nb,ndir) 
    8641085      !!--------------------------------------------- 
    8651086      !!   *** ROUTINE interpvn *** 
    8661087      !!---------------------------------------------     
    8671088      ! 
    868       INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 
    869       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 
     1089      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 
     1090      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 
    8701091      LOGICAL, INTENT(in) :: before 
     1092      INTEGER, INTENT(in) :: nb , ndir 
    8711093      ! 
    8721094      INTEGER :: ji,jj,jk 
    873       REAL(wp) :: zrhox  
     1095      REAL(wp) :: zrhox 
     1096! VERTICAL REFINEMENT BEGIN 
     1097      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child 
     1098      REAL(wp), DIMENSION(k1:k2) :: tabin 
     1099      REAL(wp) :: h_in(k1:k2) 
     1100      REAL(wp) :: h_out(1:jpk) 
     1101      INTEGER :: N_in, N_out 
     1102      REAL(wp) :: h_diff 
     1103      LOGICAL :: northern_side,southern_side 
     1104      INTEGER :: jref 
     1105 
     1106! VERTICAL REFINEMENT END 
    8741107      !!---------------------------------------------     
    8751108      !       
    … …  
    8791112            DO jj=j1,j2 
    8801113               DO ji=i1,i2 
    881                   ptab(ji,jj,jk) = e1v(ji,jj) * vn(ji,jj,jk) 
    882                   ptab(ji,jj,jk) = ptab(ji,jj,jk) * fse3v(ji,jj,jk) 
     1114                  ptab(ji,jj,jk,1) = e1v(ji,jj) * vn(ji,jj,jk) 
     1115                  ptab(ji,jj,jk,1) = ptab(ji,jj,jk,1) * fse3v(ji,jj,jk) 
     1116                  ptab(ji,jj,jk,2) = fse3v(ji,jj,jk) 
    8831117               END DO 
    8841118            END DO 
    8851119         END DO 
    886       ELSE           
     1120      ELSE         
     1121! VERTICAL REFINEMENT BEGIN 
     1122         ptab_child(:,:,:) = 0. 
     1123         southern_side = (nb == 2).AND.(ndir == 1) 
     1124         northern_side = (nb == 2).AND.(ndir == 2) 
     1125         do jj=j1,j2 
     1126         jref = jj 
     1127         IF (southern_side) jref = 2 
     1128         IF (northern_side) jref = nlcj-2 
     1129         do ji=i1,i2 
     1130 
     1131         N_in = 0 
     1132         do jk=k1,k2 
     1133           if (interp_scales_v(ji,jj,jk) == 0) EXIT 
     1134             N_in = N_in + 1 
     1135             tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 
     1136             h_in(N_in) = ptab(ji,jj,jk,2) 
     1137         enddo 
     1138         IF (N_in == 0) THEN 
     1139           ptab_child(ji,jj,:) = 0. 
     1140           CYCLE 
     1141         ENDIF 
     1142          
     1143         N_out = 0 
     1144         do jk=1,jpk 
     1145           if (vmask(ji,jref,jk) == 0) EXIT 
     1146           N_out = N_out + 1 
     1147           h_out(N_out) = fse3v(ji,jj,jk) 
     1148         enddo 
     1149         IF (N_out == 0) THEN 
     1150           ptab_child(ji,jj,:) = 0. 
     1151           CYCLE 
     1152         ENDIF 
     1153          
     1154         h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     1155         IF (h_diff > 0.) THEN 
     1156           N_in = N_in + 1 
     1157           h_in(N_in) = h_diff 
     1158           tabin(N_in) = 0. 
     1159         ELSE 
     1160           h_out(N_out) = h_out(N_out) - h_diff 
     1161         ENDIF 
     1162          
     1163         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) 
     1164          
     1165         ptab_child(ji,jj,N_out) = ptab_child(ji,jj,N_out) * h_out(N_out) / fse3v(ji,jj,N_out) 
     1166 
     1167         enddo 
     1168         enddo 
     1169! in the following 
     1170! remove division of va by fs e3v (already done) 
     1171! VERTICAL REFINEMENT END 
    8871172         zrhox= Agrif_Rhox() 
    8881173         DO jk=1,jpkm1 
    8891174            DO jj=j1,j2 
    890                va(i1:i2,jj,jk) = (ptab(i1:i2,jj,jk)/(zrhox*e1v(i1:i2,jj))) 
    891                va(i1:i2,jj,jk) = va(i1:i2,jj,jk) / fse3v(i1:i2,jj,jk) 
     1175               va(i1:i2,jj,jk) = (ptab_child(i1:i2,jj,jk)/(zrhox*e1v(i1:i2,jj))) 
    8921176            END DO 
    8931177         END DO 

  • branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90

    r5819 r6258  
    1 #define TWO_WAY        /* TWO WAY NESTING */ 
     1#undef TWO_WAY        /* TWO WAY NESTING */ 
    22#undef DECAL_FEEDBACK /* SEPARATION of INTERFACES*/ 
    33  
    … …  
    1818 
    1919   PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn,Update_Scales 
     20    
     21! VERTICAL REFINEMENT BEGIN    
     22   PUBLIC Agrif_Init_UpdateScales 
     23   REAL, DIMENSION(:,:,:), ALLOCATABLE :: update_scales_t, update_scales_u, update_scales_v 
     24!$AGRIF_DO_NOT_TREAT 
     25   LOGICAL :: scaleT, scaleU, scaleV = .FALSE. 
     26!$AGRIF_END_DO_NOT_TREAT 
     27! VERTICAL REFINEMENT END 
     28 
    2029# if defined key_zdftke 
    2130   PUBLIC Agrif_Update_Tke 
    … …  
    2332   !!---------------------------------------------------------------------- 
    2433   !! NEMO/NST 3.6 , NEMO Consortium (2010) 
    25    !! $Id$ 
     34   !! $Id: agrif_opa_update.F90 5081 2015-02-13 09:51:27Z smasson $ 
    2635   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    2736   !!---------------------------------------------------------------------- 
    2837 
     38#  include "domzgr_substitute.h90"   
     39 
    2940CONTAINS 
     41 
     42! VERTICAL REFINEMENT BEGIN 
     43   SUBROUTINE Agrif_Init_UpdateScales 
     44 
     45    scaleT = .TRUE. 
     46    Agrif_UseSpecialValueInUpdate = .FALSE. 
     47    Agrif_SpecialValueFineGrid = 0. 
     48    Call Agrif_Update_Variable(scales_t_id,procname=updatescales) 
     49    Agrif_UseSpecialValueInUpdate = .FALSE. 
     50    scaleT = .FALSE. 
     51 
     52    scaleU = .TRUE. 
     53    Call Agrif_Update_Variable(scales_u_id,procname=updatescales) 
     54    scaleU = .FALSE. 
     55 
     56    scaleV = .TRUE. 
     57    Call Agrif_Update_Variable(scales_v_id,procname=updatescales) 
     58    scaleV = .FALSE. 
     59 
     60   END SUBROUTINE Agrif_Init_UpdateScales 
     61    
     62   SUBROUTINE updatescales( ptab, i1, i2, j1, j2, k1, k2, before ) 
     63 
     64      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 
     65      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 
     66      LOGICAL, iNTENT(in) :: before 
     67 
     68      INTEGER :: ji,jj,jk 
     69 
     70      IF (before) THEN 
     71        IF (scaleT) THEN 
     72         DO jk=k1,k2 
     73            DO jj=j1,j2 
     74               DO ji=i1,i2 
     75                  ptab(ji,jj,jk) = fse3t(ji,jj,jk)*tmask(ji,jj,jk) 
     76               END DO 
     77            END DO 
     78         END DO 
     79        ELSEIF (scaleU) THEN 
     80        DO jk=k1,k2 
     81            DO jj=j1,j2 
     82               DO ji=i1,i2 
     83                  ptab(ji,jj,jk) = fse3u(ji,jj,jk)*umask(ji,jj,jk) 
     84               END DO 
     85            END DO 
     86         END DO 
     87        ELSEIF (scaleV) THEN 
     88        DO jk=k1,k2 
     89            DO jj=j1,j2 
     90               DO ji=i1,i2 
     91                  ptab(ji,jj,jk) = fse3v(ji,jj,jk)*vmask(ji,jj,jk) 
     92               END DO 
     93            END DO 
     94         END DO 
     95        ENDIF 
     96      ELSE 
     97         IF (scaleT) THEN 
     98           Allocate(update_scales_t(i1:i2,j1:j2,k1:k2)) 
     99           update_scales_t(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2) 
     100         ELSEIF (scaleU) THEN 
     101           Allocate(update_scales_u(i1:i2,j1:j2,k1:k2)) 
     102           update_scales_u(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2) 
     103         ELSEIF (scaleV) THEN 
     104           Allocate(update_scales_v(i1:i2,j1:j2,k1:k2)) 
     105           update_scales_v(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2) 
     106         ENDIF 
     107      ENDIF 
     108 
     109   END SUBROUTINE updatescales 
     110! VERTICAL REFINEMENT END 
    30111 
    31112   RECURSIVE SUBROUTINE Agrif_Update_Tra( ) 
    … …  
    160241      INTEGER, INTENT(in) :: kt 
    161242      !        
     243      return 
     244       
    162245      IF( ( Agrif_NbStepint() .NE. 0 ) .AND. (kt /= 0) ) RETURN 
    163246#  if defined TWO_WAY 
    … …  
    177260# endif /* key_zdftke */ 
    178261 
    179    SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
     262   SUBROUTINE updateTS( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    180263      !!--------------------------------------------- 
    181264      !!           *** ROUTINE updateT *** 
    182265      !!--------------------------------------------- 
    183 #  include "domzgr_substitute.h90" 
     266 
    184267      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 
    185       REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     268      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: ptab 
    186269      LOGICAL, INTENT(in) :: before 
    187270      !! 
    188271      INTEGER :: ji,jj,jk,jn 
     272! VERTICAL REFINEMENT BEGIN 
     273      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child 
     274      REAL(wp) :: h_in(k1:k2) 
     275      REAL(wp) :: h_out(1:jpk) 
     276      INTEGER :: N_in, N_out 
     277      REAL(wp) :: h_diff 
     278      REAL(wp) :: tabin(k1:k2,n1:n2) 
     279! VERTICAL REFINEMENT END 
    189280      !!--------------------------------------------- 
    190281      ! 
    … …  
    194285               DO jj=j1,j2 
    195286                  DO ji=i1,i2 
    196                      tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) 
     287                     ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) 
    197288                  END DO 
    198289               END DO 
    … …  
    200291         END DO 
    201292      ELSE 
     293! VERTICAL REFINEMENT BEGIN 
     294         ptab_child(:,:,:,:) = 0. 
     295          
     296         DO jj=j1,j2 
     297         DO ji=i1,i2 
     298           N_in = 0 
     299           DO jk=k1,k2 
     300             IF (update_scales_t(ji,jj,jk) == 0) EXIT 
     301             N_in = N_in + 1 
     302             tabin(jk,:) = ptab(ji,jj,jk,:) 
     303             h_in(N_in) = update_scales_t(ji,jj,jk) 
     304           ENDDO 
     305           N_out = 0 
     306           DO jk=1,jpk 
     307             IF (tmask(ji,jj,jk) == 0) EXIT 
     308             N_out = N_out + 1 
     309             h_out(N_out) = fse3t(ji,jj,jk) 
     310           ENDDO 
     311           IF (N_in > 0) THEN 
     312             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     313             IF (h_diff > 0) THEN 
     314               N_in = N_in + 1 
     315               h_in(N_in) = h_diff 
     316               tabin(N_in,:) = tabin(N_in-1,:) 
     317             ELSEIF (h_diff < 0) THEN 
     318             print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out)) 
     319             print *,'Nval = ',N_out,mbathy(ji,jj) 
     320             print *,'BATHY = ',gdepw_0(ji,jj,mbathy(ji,jj)+1),sum(e3t_0(ji,jj,1:mbathy(ji,jj))) 
     321          !   STOP 
     322               N_out = N_out + 1 
     323               h_out(N_out) = - h_diff 
     324             ENDIF 
     325             DO jn=n1,n2 
     326               CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),ptab_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 
     327             ENDDO 
     328          ENDIF 
     329         ENDDO 
     330         ENDDO 
     331          
     332! WARNING : 
     333! ptab replaced by ptab_child in the following 
     334! k1 k2 replaced by 1 jpk 
     335! VERTICAL REFINEMENT END 
     336 
    202337         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    203338            ! Add asselin part 
    204339            DO jn = n1,n2 
    205                DO jk=k1,k2 
     340               DO jk=1,jpk 
    206341                  DO jj=j1,j2 
    207342                     DO ji=i1,i2 
    208                         IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 
     343                        IF( ptab_child(ji,jj,jk,jn) .NE. 0. ) THEN 
    209344                           tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) &  
    210                                  & + atfp * ( tabres(ji,jj,jk,jn) & 
     345                                 & + atfp * ( ptab_child(ji,jj,jk,jn) & 
    211346                                 &             - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
    212347                        ENDIF 
    … …  
    217352         ENDIF 
    218353         DO jn = n1,n2 
    219             DO jk=k1,k2 
     354            DO jk=1,jpk 
    220355               DO jj=j1,j2 
    221356                  DO ji=i1,i2 
    222                      IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN  
    223                         tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) * tmask(ji,jj,jk) 
     357                     IF( ptab_child(ji,jj,jk,jn) .NE. 0. ) THEN  
     358                        tsn(ji,jj,jk,jn) = ptab_child(ji,jj,jk,jn) * tmask(ji,jj,jk) 
    224359                     END IF 
    225360                  END DO 
    … …  
    231366   END SUBROUTINE updateTS 
    232367 
    233    SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, before ) 
     368   SUBROUTINE updateu( ptab, i1, i2, j1, j2, k1, k2, before ) 
    234369      !!--------------------------------------------- 
    235370      !!           *** ROUTINE updateu *** 
    236371      !!--------------------------------------------- 
    237 #  include "domzgr_substitute.h90" 
    238372      !! 
    239373      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2 
    240       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres 
     374      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 
    241375      LOGICAL, INTENT(in) :: before 
    242376      !!  
    243377      INTEGER :: ji, jj, jk 
    244378      REAL(wp) :: zrhoy 
     379! VERTICAL REFINEMENT BEGIN 
     380      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child 
     381      REAL(wp) :: h_in(k1:k2) 
     382      REAL(wp) :: h_out(1:jpk) 
     383      INTEGER :: N_in, N_out 
     384      REAL(wp) :: h_diff 
     385      REAL(wp) :: tabin(k1:k2) 
     386! VERTICAL REFINEMENT END 
    245387      !!--------------------------------------------- 
    246388      !  
    … …  
    250392            DO jj=j1,j2 
    251393               DO ji=i1,i2 
    252                   tabres(ji,jj,jk) = e2u(ji,jj) * un(ji,jj,jk) 
    253                   tabres(ji,jj,jk) = tabres(ji,jj,jk) * fse3u_n(ji,jj,jk) 
    254                END DO 
    255             END DO 
    256          END DO 
    257          tabres = zrhoy * tabres 
    258       ELSE 
    259          DO jk=k1,k2 
    260             DO jj=j1,j2 
    261                DO ji=i1,i2 
    262                   tabres(ji,jj,jk) = tabres(ji,jj,jk) / e2u(ji,jj) / fse3u_n(ji,jj,jk) 
     394                  ptab(ji,jj,jk) = e2u(ji,jj) * un(ji,jj,jk) 
     395                  ptab(ji,jj,jk) = ptab(ji,jj,jk) * fse3u_n(ji,jj,jk) 
     396               END DO 
     397            END DO 
     398         END DO 
     399         ptab = zrhoy * ptab 
     400      ELSE 
     401! VERTICAL REFINEMENT BEGIN 
     402         ptab_child(:,:,:) = 0. 
     403          
     404         DO jj=j1,j2 
     405         DO ji=i1,i2 
     406           N_in = 0 
     407           DO jk=k1,k2 
     408             IF (update_scales_u(ji,jj,jk) == 0) EXIT 
     409             N_in = N_in + 1 
     410             tabin(jk) = ptab(ji,jj,jk)/update_scales_u(ji,jj,jk) 
     411             h_in(N_in) = update_scales_u(ji,jj,jk) 
     412           ENDDO 
     413           N_out = 0 
     414           DO jk=1,jpk 
     415             IF (umask(ji,jj,jk) == 0) EXIT 
     416             N_out = N_out + 1 
     417             h_out(N_out) = fse3u(ji,jj,jk) 
     418           ENDDO 
     419           IF (N_in * N_out > 0) THEN 
     420             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     421             if (h_diff < 0.) then 
     422             print *,'CHECK YOUR BATHY ...' 
     423             stop 
     424             else ! Extends with 0 
     425             N_in = N_in + 1 
     426             tabin(N_in) = 0. 
     427             h_in(N_in) = h_diff 
     428             endif 
     429             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) 
     430          ENDIF 
     431         ENDDO 
     432         ENDDO 
     433          
     434! WARNING : 
     435! ptab replaced by ptab_child in the following 
     436! k1 k2 replaced by 1 jpk 
     437! remove division by fs e3u (already done) 
     438! VERTICAL REFINEMENT END 
     439 
     440         DO jk=1,jpk 
     441            DO jj=j1,j2 
     442               DO ji=i1,i2 
     443                  ptab_child(ji,jj,jk) = ptab_child(ji,jj,jk) / e2u(ji,jj) 
    263444                  ! 
    264445                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    265446                     ub(ji,jj,jk) = ub(ji,jj,jk) &  
    266                            & + atfp * ( tabres(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 
     447                           & + atfp * ( ptab_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 
    267448                  ENDIF 
    268449                  ! 
    269                   un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk) 
     450                  un(ji,jj,jk) = ptab_child(ji,jj,jk) * umask(ji,jj,jk) 
    270451               END DO 
    271452            END DO 
    … …  
    275456   END SUBROUTINE updateu 
    276457 
    277    SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, before ) 
     458   SUBROUTINE updatev( ptab, i1, i2, j1, j2, k1, k2, before ) 
    278459      !!--------------------------------------------- 
    279460      !!           *** ROUTINE updatev *** 
    280461      !!--------------------------------------------- 
    281 #  include "domzgr_substitute.h90" 
    282462      !! 
    283463      INTEGER :: i1,i2,j1,j2,k1,k2 
    284464      INTEGER :: ji,jj,jk 
    285       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: tabres 
     465      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab 
    286466      LOGICAL :: before 
    287467      !! 
    288468      REAL(wp) :: zrhox 
     469! VERTICAL REFINEMENT BEGIN 
     470      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child 
     471      REAL(wp) :: h_in(k1:k2) 
     472      REAL(wp) :: h_out(1:jpk) 
     473      INTEGER :: N_in, N_out 
     474      REAL(wp) :: h_diff 
     475      REAL(wp) :: tabin(k1:k2) 
     476! VERTICAL REFINEMENT END 
    289477      !!---------------------------------------------       
    290478      ! 
    … …  
    294482            DO jj=j1,j2 
    295483               DO ji=i1,i2 
    296                   tabres(ji,jj,jk) = e1v(ji,jj) * vn(ji,jj,jk) 
    297                   tabres(ji,jj,jk) = tabres(ji,jj,jk) * fse3v_n(ji,jj,jk) 
    298                END DO 
    299             END DO 
    300          END DO 
    301          tabres = zrhox * tabres 
    302       ELSE 
     484                  ptab(ji,jj,jk) = e1v(ji,jj) * vn(ji,jj,jk) 
     485                  ptab(ji,jj,jk) = ptab(ji,jj,jk) * fse3v_n(ji,jj,jk) 
     486               END DO 
     487            END DO 
     488         END DO 
     489         ptab = zrhox * ptab 
     490      ELSE 
     491! VERTICAL REFINEMENT BEGIN 
     492         ptab_child(:,:,:) = 0. 
     493          
     494         DO jj=j1,j2 
     495         DO ji=i1,i2 
     496           N_in = 0 
     497           DO jk=k1,k2 
     498             IF (update_scales_v(ji,jj,jk) == 0) EXIT 
     499             N_in = N_in + 1 
     500             tabin(jk) = ptab(ji,jj,jk)/update_scales_v(ji,jj,jk) 
     501             h_in(N_in) = update_scales_v(ji,jj,jk) 
     502           ENDDO 
     503           N_out = 0 
     504           DO jk=1,jpk 
     505             IF (vmask(ji,jj,jk) == 0) EXIT 
     506             N_out = N_out + 1 
     507             h_out(N_out) = fse3v(ji,jj,jk) 
     508           ENDDO 
     509           IF (N_in * N_out > 0) THEN 
     510             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     511             if (h_diff < 0.) then 
     512             print *,'CHECK YOUR BATHY ...' 
     513             stop 
     514             else ! Extends with 0 
     515             N_in = N_in + 1 
     516             tabin(N_in) = 0. 
     517             h_in(N_in) = h_diff 
     518             endif 
     519             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) 
     520          ENDIF 
     521         ENDDO 
     522         ENDDO 
     523          
     524! WARNING : 
     525! ptab replaced by ptab_child in the following 
     526! k1 k2 replaced by 1 jpk 
     527! remove division by fs e3v (already done) 
     528! VERTICAL REFINEMENT END 
     529 
    303530         DO jk=k1,k2 
    304531            DO jj=j1,j2 
    305532               DO ji=i1,i2 
    306                   tabres(ji,jj,jk) = tabres(ji,jj,jk) / e1v(ji,jj) / fse3v_n(ji,jj,jk) 
     533                  ptab_child(ji,jj,jk) = ptab_child(ji,jj,jk) / e1v(ji,jj) 
    307534                  ! 
    308535                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    309536                     vb(ji,jj,jk) = vb(ji,jj,jk) &  
    310                            & + atfp * ( tabres(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 
     537                           & + atfp * ( ptab_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 
    311538                  ENDIF 
    312539                  ! 
    313                   vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk) 
     540                  vn(ji,jj,jk) = ptab_child(ji,jj,jk) * vmask(ji,jj,jk) 
    314541               END DO 
    315542            END DO 

  • branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_user.F90

    r5819 r6258  
    22!!---------------------------------------------------------------------- 
    33!! NEMO/NST 3.4 , NEMO Consortium (2012) 
    4 !! $Id$ 
     4!! $Id: agrif_user.F90 5081 2015-02-13 09:51:27Z smasson $ 
    55!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    66!!---------------------------------------------------------------------- 
    … …  
    1717   USE par_oce 
    1818   USE dom_oce 
     19   USE Agrif_Util 
     20   USE lib_mpp         ! distributed memory computing 
    1921   USE nemogcm 
    2022   ! 
    2123   IMPLICIT NONE 
     24   INTEGER ::   ios 
     25   LOGICAL ::  is_open 
     26   NAMELIST/namcfg/ cp_cfg, cp_cfz, jp_cfg, jpidta, jpjdta, jpkdta, jpiglo, jpjglo, & 
     27         &             jpizoom, jpjzoom, jperio 
    2228   !!---------------------------------------------------------------------- 
    2329   ! 
    … …  
    4551      nperio  = 0 
    4652      jperio  = 0 
     53   ELSE 
     54   IF (Agrif_Nb_step() == 0) THEN 
     55      INQUIRE(file = 'namelist_ref', opened = is_open) 
     56      IF (.not.is_open) THEN 
     57      !                             ! Open reference namelist and configuration namelist files 
     58      CALL ctl_opn( numnam_ref, 'namelist_ref', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. ) 
     59      CALL ctl_opn( numnam_cfg, 'namelist_cfg', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE. ) 
     60       
     61      REWIND( numnam_ref )              ! Namelist namcfg in reference namelist : Control prints & Benchmark 
     62      READ  ( numnam_ref, namcfg, IOSTAT = ios, ERR = 903 ) 
     63903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfg in reference namelist', .TRUE. ) 
     64 
     65      REWIND( numnam_cfg )              ! Namelist namcfg in confguration namelist : Control prints & Benchmark 
     66      READ  ( numnam_cfg, namcfg, IOSTAT = ios, ERR = 904 ) 
     67904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcfg in configuration namelist', .TRUE. )  
     68 
     69      IF( numnam_ref      /= -1 )   CLOSE( numnam_ref      )   ! oce reference namelist 
     70      IF( numnam_cfg      /= -1 )   CLOSE( numnam_cfg      )   ! oce configuration namelist 
     71      ENDIF 
     72   ENDIF 
    4773   ENDIF 
    4874   ! 
    … …  
    6793   ! 0. Initializations 
    6894   !------------------- 
    69    IF( cp_cfg == 'orca' ) THEN 
    70       IF ( jp_cfg == 2 .OR. jp_cfg == 025 .OR. jp_cfg == 05 & 
    71             &                      .OR. jp_cfg == 4 ) THEN 
    72          jp_cfg = -1    ! set special value for jp_cfg on fine grids 
    73          cp_cfg = "default" 
    74       ENDIF 
    75    ENDIF 
    7695   ! Specific fine grid Initializations 
    7796   ! no tracer damping on fine grids 
    … …  
    188207   !--------------------------------------------------------------------- 
    189208   CALL agrif_declare_var 
     209 
     210! VERTICAL REFINEMENT BEGIN 
     211   CALL Agrif_Init_InterpScales() 
     212   CALL Agrif_Init_UpdateScales() 
     213! VERTICAL REFINEMENT END 
    190214 
    191215   ! 2. First interpolations of potentially non zero fields 
    … …  
    289313         CALL Agrif_Bc_variable(vmsk_id,calledweight=1.,procname=interpvmsk) 
    290314    ! check if tmask and vertical scale factors agree with parent over first two coarse grid points: 
    291          CALL Agrif_Bc_variable(e3t_id,calledweight=1.,procname=interpe3t) 
     315!        CALL Agrif_Bc_variable(e3t_id,calledweight=1.,procname=interpe3t) 
    292316         ! 
    293317         IF (lk_mpp) CALL mpp_sum( kindic_agr ) 
    … …  
    343367   CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_sponge_id) 
    344368 
    345    CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_interp_id) 
    346    CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_interp_id) 
     369   CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_interp_id) 
     370   CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_interp_id) 
    347371   CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_update_id) 
    348372   CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_update_id) 
    … …  
    353377   CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),umsk_id) 
    354378   CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vmsk_id) 
    355  
    356    CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,3/),scales_t_id) 
    357379 
    358380   CALL agrif_declare_variable((/1,2/),(/2,3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),unb_id) 
    … …  
    365387   CALL agrif_declare_variable((/2,2/),(/3,3/),(/'x','y'/),(/1,1/),(/nlci,nlcj/),sshn_id) 
    366388 
     389! VERTICAL REFINEMENT BEGIN 
     390   CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),scales_t_id) 
     391   CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),scales_u_id) 
     392   CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),scales_v_id) 
     393! VERTICAL REFINEMENT END 
     394 
    367395# if defined key_zdftke 
    368396   CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/), en_id) 
    … …  
    393421   CALL Agrif_Set_bcinterp(umsk_id,interp=AGRIF_constant) 
    394422   CALL Agrif_Set_bcinterp(vmsk_id,interp=AGRIF_constant) 
     423 
     424! VERTICAL REFINEMENT BEGIN 
     425   CALL Agrif_Set_bcinterp(scales_t_id,interp=Agrif_linear) 
     426   CALL Agrif_Set_bcinterp(scales_u_id,interp1=Agrif_linear,interp2=Agrif_constant) 
     427   CALL Agrif_Set_bcinterp(scales_v_id,interp1=Agrif_constant,interp2=AGRIF_linear) 
     428! VERTICAL REFINEMENT END 
    395429 
    396430# if defined key_zdftke 
    … …  
    422456   CALL Agrif_Set_bc(vmsk_id,(/0,0/)) 
    423457 
     458! VERTICAL REFINEMENT BEGIN 
     459   CALL Agrif_Set_bc(scales_t_id,(/-nn_sponge_len*Agrif_irhox()-1,1/)) 
     460   CALL Agrif_Set_bc(scales_u_id,(/-nn_sponge_len*Agrif_irhox()-1,1/)) 
     461   CALL Agrif_Set_bc(scales_v_id,(/-nn_sponge_len*Agrif_irhox()-1,1/)) 
     462! VERTICAL REFINEMENT END 
     463 
    424464# if defined key_zdftke 
    425465   CALL Agrif_Set_bc(avm_id ,(/0,1/)) 
    … …  
    439479   CALL Agrif_Set_Updatetype(ub2b_update_id,update1 = Agrif_Update_Copy, update2 = Agrif_Update_Average) 
    440480   CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy) 
     481 
     482! VERTICAL REFINEMENT BEGIN 
     483   CALL Agrif_Set_Updatetype(scales_t_id, update = AGRIF_Update_Average) 
     484   CALL Agrif_Set_Updatetype(scales_u_id,update1 = Agrif_Update_Copy, update2 = Agrif_Update_Average) 
     485   CALL Agrif_Set_Updatetype(scales_v_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy) 
     486! VERTICAL REFINEMENT END 
    441487 
    442488# if defined key_zdftke 

  • branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r5819 r6258  
    405405      !! 
    406406      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout 
    407       CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names 
     407      CHARACTER (len=60) ::   clhstnam, clop, clmx           ! local names 
    408408      INTEGER  ::   inum = 11                                ! temporary logical unit 
    409409      INTEGER  ::   ji, jj, jk                               ! dummy loop indices 
    … …  
    779779      IF( lk_vvl ) THEN 
    780780         zw3d(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 
    781          CALL histwrite( nid_T, "vovvle3t", it, fse3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness 
     781         CALL histwrite( nid_T, "vovvle3t", it, fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! level thickness 
    782782         CALL histwrite( nid_T, "vovvldep", it, fsdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth 
    783783         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation 

  • branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r5819 r6258  
    9797      END DO 
    9898      ! 
    99       IF( lk_vvl )           CALL dom_vvl_init ! Vertical variable mesh 
     99      IF( lk_vvl )           THEN 
     100        CALL dom_vvl_init ! Vertical variable mesh 
     101      ELSE 
     102        ln_vvl_ztilde = .FALSE. 
     103        ln_vvl_layer  = .FALSE. 
     104      ENDIF 
     105 
    100106      ! 
    101107      IF( lk_c1d         )   CALL cor_c1d      ! 1D configuration: Coriolis set at T-point 

  • branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r5819 r6258  
    620620      ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt 
    621621      ! - JC - hu_b, hv_b, hur_b, hvr_b also 
    622       CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F'  ) 
     622      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F'  ) 
    623623      ! Vertical scale factor interpolations 
    624624      ! ------------------------------------ 
    625       CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  ) 
     625      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n(:,:,:), 'W'  ) 
    626626      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 
    627627      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 
    628       CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W'  ) 
     628      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b(:,:,:), 'W'  ) 
    629629      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 
    630630      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 
    … …  
    652652      ! Local depth and Inverse of the local depth of the water column at u- and v- points 
    653653      ! ---------------------------------------------------------------------------------- 
    654       hu (:,:) = hu_a (:,:) 
    655       hv (:,:) = hv_a (:,:) 
     654      hu (:,:) = hu_a(:,:) 
     655      hv (:,:) = hv_a(:,:) 
    656656 
    657657      ! Inverse of the local depth 
    … …  
    668668      ! Write outputs 
    669669      ! ============= 
    670       CALL iom_put(     "e3t" , fse3t_n  (:,:,:) ) 
    671       CALL iom_put(     "e3u" , fse3u_n  (:,:,:) ) 
    672       CALL iom_put(     "e3v" , fse3v_n  (:,:,:) ) 
    673       CALL iom_put(     "e3w" , fse3w_n  (:,:,:) ) 
    674       CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) ) 
     670      CALL iom_put(     "e3t" , fse3t_n(:,:,:) ) 
     671      CALL iom_put(     "e3u" , fse3u_n(:,:,:) ) 
     672      CALL iom_put(     "e3v" , fse3v_n(:,:,:) ) 
     673      CALL iom_put(     "e3w" , fse3w_n(:,:,:) ) 
     674      CALL iom_put( "tpt_dep" , fsde3w_n(:,:,:) ) 
    675675      IF( iom_use("e3tdef") )   & 
    676676         CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 

  • branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5819 r6258  
    294294               e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   ) 
    295295            END DO 
     296! MODIFY THE COMPUTATION OF GDEPW_1D - Laurent Debreu Change for Agrif Vertical interpolation 
     297            gdepw_1d(1) = 0. 
     298            do jk=2,jpk 
     299            gdepw_1d(jk) = gdepw_1d(jk-1)+e3t_1d(jk-1) 
     300            enddo 
    296301         ELSE 
    297302            DO jk = 1, jpk 
    … …  
    898903      ! 
    899904      DO jk = 1, jpk 
    900          gdept_0 (:,:,jk) = gdept_1d(jk) 
    901          gdepw_0 (:,:,jk) = gdepw_1d(jk) 
     905         gdept_0(:,:,jk) = gdept_1d(jk) 
     906         gdepw_0(:,:,jk) = gdepw_1d(jk) 
    902907         gdep3w_0(:,:,jk) = gdepw_1d(jk) 
    903          e3t_0   (:,:,jk) = e3t_1d  (jk) 
    904          e3u_0   (:,:,jk) = e3t_1d  (jk) 
    905          e3v_0   (:,:,jk) = e3t_1d  (jk) 
    906          e3f_0   (:,:,jk) = e3t_1d  (jk) 
    907          e3w_0   (:,:,jk) = e3w_1d  (jk) 
    908          e3uw_0  (:,:,jk) = e3w_1d  (jk) 
    909          e3vw_0  (:,:,jk) = e3w_1d  (jk) 
     908         e3t_0(:,:,jk) = e3t_1d(jk) 
     909         e3u_0(:,:,jk) = e3t_1d(jk) 
     910         e3v_0(:,:,jk) = e3t_1d(jk) 
     911         e3f_0(:,:,jk) = e3t_1d(jk) 
     912         e3w_0(:,:,jk) = e3w_1d(jk) 
     913         e3uw_0(:,:,jk) = e3w_1d(jk) 
     914         e3vw_0(:,:,jk) = e3w_1d(jk) 
    910915      END DO 
    911916      ! 
    … …  
    10431048                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
    10441049                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
    1045                   e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
     1050                  e3t_0(ji,jj,ik) = e3t_1d(ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
    10461051                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
    10471052                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
    10481053                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
    10491054                  !       ... on ik+1 
    1050                   e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1051                   e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     1055                  e3w_0(ji,jj,ik+1) = e3t_0(ji,jj,ik) 
     1056                  e3t_0(ji,jj,ik+1) = e3t_0(ji,jj,ik) 
    10521057                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
    10531058               ENDIF 
    … …  
    10891094                     &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
    10901095                     &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
    1091                   e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
    1092                   e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
     1096                  e3t_0(ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
     1097                  e3w_0(ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
    10931098 
    10941099                  IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
    1095                      e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
     1100                     e3w_0(ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
    10961101                  ENDIF  
    10971102               !       ... on ik / ik-1  
    … …  
    21212126      fsdepw(:,:,:) = gdepw_0 (:,:,:) 
    21222127      fsde3w(:,:,:) = gdep3w_0(:,:,:) 
    2123       fse3t (:,:,:) = e3t_0   (:,:,:) 
    2124       fse3u (:,:,:) = e3u_0   (:,:,:) 
    2125       fse3v (:,:,:) = e3v_0   (:,:,:) 
    2126       fse3f (:,:,:) = e3f_0   (:,:,:) 
    2127       fse3w (:,:,:) = e3w_0   (:,:,:) 
    2128       fse3uw(:,:,:) = e3uw_0  (:,:,:) 
    2129       fse3vw(:,:,:) = e3vw_0  (:,:,:) 
     2128      fse3t(:,:,:) = e3t_0(:,:,:) 
     2129      fse3u(:,:,:) = e3u_0(:,:,:) 
     2130      fse3v(:,:,:) = e3v_0(:,:,:) 
     2131      fse3f(:,:,:) = e3f_0(:,:,:) 
     2132      fse3w(:,:,:) = e3w_0(:,:,:) 
     2133      fse3uw(:,:,:) = e3uw_0(:,:,:) 
     2134      fse3vw(:,:,:) = e3vw_0(:,:,:) 
    21302135!! 
    21312136      ! HYBRID :  

  • branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r5819 r6258  
    845845               &                   * (  rhd(ji,jj,1)                                    & 
    846846               &                     + 0.5_wp * ( rhd(ji,jj,2) - rhd(ji,jj,1) )         & 
    847                &                              * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) )   & 
     847               &                              * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) )   & 
    848848               &                              / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) )  ) 
    849849         END DO 

  • branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90

    r5819 r6258  
    127127         ! read filtered free surface arrays in restart file 
    128128         ! when using agrif, sshn, gcx have to be read in istate 
    129          IF(.NOT. lk_agrif)   CALL flt_rst( nit000, 'READ' )      ! read or initialize the following fields: 
     129         IF(.NOT. lk_agrif) THEN 
     130            CALL flt_rst( nit000, 'READ' )      ! read or initialize the following fields: 
    130131         !                                                        ! gcx, gcxb 
     132         ELSE 
     133            gcx (:,:) = 0.e0 
     134            gcxb(:,:) = 0.e0 
     135         ENDIF 
    131136      ENDIF 
    132137 

  • branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/SOL/solmat.F90

    r5819 r6258  
    213213      DO jj = 1, jpj 
    214214         DO ji = 1, jpi 
    215             IF( bmask(ji,jj) /= 0.e0 )   gcdprc(ji,jj) = 1.e0 / gcdmat(ji,jj) 
     215            IF( bmask(ji,jj) /= 0.e0 .AND.(gcdmat(ji,jj)/=0.) )   gcdprc(ji,jj) = 1.e0 / gcdmat(ji,jj) 
    216216         END DO 
    217217      END DO 

  • branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90

    r5819 r6258  
    210210               DO ji = fs_2, fs_jpim1 
    211211                  ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,jk) 
    212                   ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t  (ji,jj,jk) 
     212                  ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,jk) 
    213213                  zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn)   ! zrhs=right hand side  
    214214                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn) 

  • branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/TRD/trdglo.F90

    r5819 r6258  
    8585      IF( MOD(kt,nn_trd) == 0 .OR. kt == nit000 .OR. kt == nitend ) THEN 
    8686         ! 
    87          SELECT CASE( ctype ) 
    88          ! 
    89          CASE( 'TRA' )          !==  Tracers (T & S)  ==! 
    90             DO jk = 1, jpkm1       ! global sum of mask volume trend and trend*T (including interior mask) 
    91                DO jj = 1, jpj 
    92                   DO ji = 1, jpi         
    93                      zvm = e1e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 
    94                      zvt = ptrdx(ji,jj,jk) * zvm 
    95                      zvs = ptrdy(ji,jj,jk) * zvm 
    96                      tmo(ktrd) = tmo(ktrd) + zvt    
    97                      smo(ktrd) = smo(ktrd) + zvs 
    98                      t2 (ktrd) = t2(ktrd)  + zvt * tsn(ji,jj,jk,jp_tem) 
    99                      s2 (ktrd) = s2(ktrd)  + zvs * tsn(ji,jj,jk,jp_sal) 
     87      SELECT CASE( ctype ) 
     88      ! 
     89      CASE( 'TRA' )          !==  Tracers (T & S)  ==! 
     90         DO jk = 1, jpkm1       ! global sum of mask volume trend and trend*T (including interior mask) 
     91            DO jj = 1, jpj 
     92               DO ji = 1, jpi         
     93                  zvm = e1e2t(ji,jj) * fse3t(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 
     94                  zvt = ptrdx(ji,jj,jk) * zvm 
     95                  zvs = ptrdy(ji,jj,jk) * zvm 
     96                  tmo(ktrd) = tmo(ktrd) + zvt    
     97                  smo(ktrd) = smo(ktrd) + zvs 
     98                  t2 (ktrd) = t2(ktrd)  + zvt * tsn(ji,jj,jk,jp_tem) 
     99                  s2 (ktrd) = s2(ktrd)  + zvs * tsn(ji,jj,jk,jp_sal) 
    100100                  END DO 
    101101               END DO 
    … …  
    113113               CALL glo_tra_wri( kt )             ! print the results in ocean.output 
    114114               !                 
    115                tmo(:) = 0._wp                     ! prepare the next time step (domain averaged array reset to zero) 
    116                smo(:) = 0._wp 
     115            tmo(:) = 0._wp                     ! prepare the next time step (domain averaged array reset to zero) 
     116            smo(:) = 0._wp 
    117117               t2 (:) = 0._wp 
    118118               s2 (:) = 0._wp 
    … …  
    348348            WRITE (numout,9541) ( hke(jpdyn_keg) + hke(jpdyn_rvo) + hke(jpdyn_zad) ) / tvolt 
    349349            WRITE (numout,9542) ( hke(jpdyn_keg) + hke(jpdyn_zad) ) / tvolt 
    350             WRITE (numout,9543) ( hke(jpdyn_pvo) ) / tvolt 
    351             WRITE (numout,9544) ( hke(jpdyn_rvo) ) / tvolt 
    352             WRITE (numout,9545) ( hke(jpdyn_spg) ) / tvolt 
    353             WRITE (numout,9546) ( hke(jpdyn_ldf) ) / tvolt 
    354             WRITE (numout,9547) ( hke(jpdyn_zdf) ) / tvolt 
    355             WRITE (numout,9548) ( hke(jpdyn_hpg) ) / tvolt, rpktrd / tvolt 
     350            WRITE (numout,9543)   hke(jpdyn_pvo)  / tvolt 
     351            WRITE (numout,9544)   hke(jpdyn_rvo)  / tvolt 
     352            WRITE (numout,9545)   hke(jpdyn_spg)  / tvolt 
     353            WRITE (numout,9546)   hke(jpdyn_ldf)  / tvolt 
     354            WRITE (numout,9547)   hke(jpdyn_zdf)  / tvolt 
     355            WRITE (numout,9548)   hke(jpdyn_hpg)  / tvolt, rpktrd / tvolt 
    356356            WRITE (numout,*) 
    357357            WRITE (numout,*) 

  • branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl.F90

    r5819 r6258  
    464464      itmod = kt - nit000 + 1 
    465465 
    466       MODULO_NTRD : IF( MOD( itmod, nn_trd ) == 0 ) THEN        ! nitend MUST be multiple of nn_trd 
     466      IF( MOD( itmod, nn_trd ) == 0 ) THEN        ! nitend MUST be multiple of nn_trd 
    467467         ! 
    468468         ztmltot (:,:) = 0.e0   ;   zsmltot (:,:) = 0.e0   ! reset arrays to zero 
    … …  
    636636         END IF 
    637637         ! 
    638       END IF MODULO_NTRD 
     638      END IF 
    639639 
    640640      ! ====================================================================== 

  • branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r5819 r6258  
    121121#if defined key_agrif 
    122122      CALL Agrif_Init_Grids()      ! AGRIF: set the meshes 
     123      nstop = 0 ! This is getting set to 2 somewhere Agrif_Init_Grids so reset to zero for now 
    123124#endif 
    124125 
Note: See TracChangeset for help on using the changeset viewer.