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/NST_SRC/agrif_opa_interp.F90 – 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.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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 
Note: See TracChangeset for help on using the changeset viewer.