Changeset 8893
- Timestamp:
- 2017-12-04T16:30:52+01:00 (7 years ago)
- Location:
- branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC
- Files:
-
- 5 edited
-
agrif_opa_sponge.F90 (modified) (3 diffs)
-
agrif_top_interp.F90 (modified) (5 diffs)
-
agrif_top_sponge.F90 (modified) (1 diff)
-
agrif_top_update.F90 (modified) (3 diffs)
-
agrif_user.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90
r8835 r8893 216 216 DO jj=j1,j2 217 217 DO ji=i1,i2 218 tabres(ji,jj,jk,jn) = ts b(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk)218 tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 219 219 END DO 220 220 END DO … … 244 244 N_out = 0 245 245 DO jk=1,jpk ! jpk of child grid 246 IF (tmask(ji,jj,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 model246 IF (tmask(ji,jj,jk) == 0) EXIT 247 247 N_out = N_out + 1 248 248 h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above … … 264 264 do jj=j1,j2 265 265 do ji=i1,i2 266 ! This will be slow - Is there a way to speed it up and avoid divide by zero? 267 IF (tabres(ji,jj,jk,n2) .NE. 0) THEN 268 tsbdiff(ji,jj,jk,n1:n2-1) = tsb(ji,jj,jk,n1:n2) - tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 269 ELSE 270 tsbdiff(ji,jj,jk,n1:n2-1) = 0._wp 271 ENDIF 266 tsbdiff(ji,jj,jk,n1:n2-1) = tsb(ji,jj,jk,n1:n2) - tmask(ji,jj,jk)*tabres(ji,jj,jk,n1:n2-1)/e3t_n(ji,jj,jk) 272 267 enddo 273 268 enddo -
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_top_interp.F90
r6140 r8893 49 49 INTEGER, INTENT(in) :: nb , ndir 50 50 ! 51 INTEGER :: ji, jj, jk, jn ! dummy loop indices52 INTEGER :: imin, imax, jmin, jmax51 INTEGER :: ji, jj, jk, jn, iref, jref ! dummy loop indices 52 INTEGER :: imin, imax, jmin, jmax, N_in, N_out 53 53 REAL(wp) :: zrhox , zalpha1, zalpha2, zalpha3 54 54 REAL(wp) :: zalpha4, zalpha5, zalpha6, zalpha7 55 55 LOGICAL :: western_side, eastern_side,northern_side,southern_side 56 56 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child 57 REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 58 REAL(wp) :: h_in(k1:k2) 59 REAL(wp) :: h_out(1:jpk) 60 REAL(wp) :: h_diff, zrhoxy 61 62 zrhoxy = Agrif_rhox()*Agrif_rhoy() 57 63 IF (before) THEN 58 ptab(i1:i2,j1:j2,k1:k2,n1:n2) = trn(i1:i2,j1:j2,k1:k2,n1:n2) 59 ELSE 60 ! 64 DO jn = n1,n2-1 65 DO jk=k1,k2 66 DO jj=j1,j2 67 DO ji=i1,i2 68 ptab(ji,jj,jk,jn) = trn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 69 END DO 70 END DO 71 END DO 72 END DO 73 DO jk=k1,k2 74 DO jj=j1,j2 75 DO ji=i1,i2 76 ptab(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 77 END DO 78 END DO 79 END DO 80 ELSE 61 81 western_side = (nb == 1).AND.(ndir == 1) 62 82 eastern_side = (nb == 1).AND.(ndir == 2) 63 83 southern_side = (nb == 2).AND.(ndir == 1) 64 84 northern_side = (nb == 2).AND.(ndir == 2) 85 #ifdef key_vertical 86 ptab_child(:,:,:,:) = 0. 87 DO jj=j1,j2 88 DO ji=i1,i2 89 iref = ji 90 jref = jj 91 if(western_side) iref=2 92 if(eastern_side) iref=nlci-1 93 if(southern_side) jref=2 94 if(northern_side) jref=nlcj-1 95 N_in = 0 96 DO jk=k1,k2 !k2 = jpk of parent grid 97 IF (ptab(ji,jj,jk,n2) == 0) EXIT 98 N_in = N_in + 1 99 tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)/(ptab(ji,jj,jk,n2)) 100 h_in(N_in) = ptab(ji,jj,jk,n2)/(e1e2t(ji,jj)*zrhoxy) 101 END DO 102 N_out = 0 103 DO jk=1,jpk ! jpk of child grid 104 IF (tmask(iref,jref,jk) == 0) EXIT 105 N_out = N_out + 1 106 h_out(jk) = e3t_n(iref,jref,jk) 107 ENDDO 108 IF (N_in > 0) THEN 109 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 110 DO jn=1,jpts 111 call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 112 ENDDO 113 ENDIF 114 ENDDO 115 ENDDO 116 #else 117 DO jk=k1,k2 118 DO jj=j1,j2 119 DO ji=i1,i2 120 ptab_child(ji,jj,jk,n1:n2-1) = ptab(ji,jj,jk,n1:n2-1)/(zrhoxy*e1e2t(ji,jj)*e3t_n(ji,jj,jk) 121 ENDDO 122 ENDDO 123 ENDDO 124 #endif 65 125 ! 66 126 zrhox = Agrif_Rhox() … … 89 149 IF( eastern_side) THEN 90 150 DO jn = 1, jptra 91 tra(nlci,j1:j2,k1:k2,jn) = zalpha1 * ptab (nlci,j1:j2,k1:k2,jn) + zalpha2 * ptab(nlci-1,j1:j2,k1:k2,jn)151 tra(nlci,j1:j2,k1:k2,jn) = zalpha1 * ptab_child(nlci,j1:j2,k1:k2,jn) + zalpha2 * ptab_child(nlci-1,j1:j2,k1:k2,jn) 92 152 DO jk = 1, jpkm1 93 153 DO jj = jmin,jmax … … 108 168 IF( northern_side ) THEN 109 169 DO jn = 1, jptra 110 tra(i1:i2,nlcj,k1:k2,jn) = zalpha1 * ptab (i1:i2,nlcj,k1:k2,jn) + zalpha2 * ptab(i1:i2,nlcj-1,k1:k2,jn)170 tra(i1:i2,nlcj,k1:k2,jn) = zalpha1 * ptab_child(i1:i2,nlcj,k1:k2,jn) + zalpha2 * ptab_child(i1:i2,nlcj-1,k1:k2,jn) 111 171 DO jk = 1, jpkm1 112 172 DO ji = imin,imax … … 127 187 IF( western_side) THEN 128 188 DO jn = 1, jptra 129 tra(1,j1:j2,k1:k2,jn) = zalpha1 * ptab (1,j1:j2,k1:k2,jn) + zalpha2 * ptab(2,j1:j2,k1:k2,jn)189 tra(1,j1:j2,k1:k2,jn) = zalpha1 * ptab_child(1,j1:j2,k1:k2,jn) + zalpha2 * ptab_child(2,j1:j2,k1:k2,jn) 130 190 DO jk = 1, jpkm1 131 191 DO jj = jmin,jmax … … 145 205 IF( southern_side ) THEN 146 206 DO jn = 1, jptra 147 tra(i1:i2,1,k1:k2,jn) = zalpha1 * ptab (i1:i2,1,k1:k2,jn) + zalpha2 * ptab(i1:i2,2,k1:k2,jn)207 tra(i1:i2,1,k1:k2,jn) = zalpha1 * ptab_child(i1:i2,1,k1:k2,jn) + zalpha2 * ptab_child(i1:i2,2,k1:k2,jn) 148 208 DO jk=1,jpk 149 209 DO ji=imin,imax -
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_top_sponge.F90
r6140 r8893 70 70 REAL(wp), DIMENSION(i1:i2,j1:j2) :: ztu, ztv 71 71 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) :: trbdiff 72 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child 73 REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 74 REAL(wp) :: h_in(k1:k2) 75 REAL(wp) :: h_out(1:jpk) 76 INTEGER :: N_in, N_out 77 REAL(wp) :: h_diff, zrhoxy 72 78 !!---------------------------------------------------------------------- 73 79 ! 74 80 IF( before ) THEN 75 tabres(i1:i2,j1:j2,k1:k2,n1:n2) = trn(i1:i2,j1:j2,k1:k2,n1:n2) 76 ELSE 77 !!gm line below use of :,: versus i1:i2,j1:j2 .... strange, not wrong. ===>> to be corrected 78 trbdiff(:,:,:,:) = trb(i1:i2,j1:j2,:,:) - tabres(:,:,:,:) 81 DO jn = n1,n2-1 82 DO jk=k1,k2 83 DO jj=j1,j2 84 DO ji=i1,i2 85 tabres(ji,jj,jk,jn) = trn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 86 END DO 87 END DO 88 END DO 89 END DO 90 DO jk=k1,k2 91 DO jj=j1,j2 92 DO ji=i1,i2 93 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 94 END DO 95 END DO 96 END DO 97 ELSE 98 ! 99 #ifdef key_vertical 100 tabres_child(:,:,:,:) = 0. 101 zrhoxy = Agrif_rhox()*Agrif_rhoy() 102 DO jj=j1,j2 103 DO ji=i1,i2 104 N_in = 0 105 DO jk=k1,k2 !k2 = jpk of parent grid 106 IF (tabres(ji,jj,jk,n2) == 0) EXIT 107 N_in = N_in + 1 108 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 109 h_in(N_in) = tabres(ji,jj,jk,n2)/(e1e2t(ji,jj)*zrhoxy) 110 END DO 111 N_out = 0 112 DO jk=1,jpk ! jpk of child grid 113 IF (tmask(ji,jj,jk) == 0) EXIT 114 N_out = N_out + 1 115 h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 116 ENDDO 117 IF (N_in > 0) THEN 118 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 119 tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for????? 120 DO jn=1,jpts 121 call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 122 ENDDO 123 ENDIF 124 DO jk=1,jpk 125 tsbdiff(ji,jj,jk,n1:n2-1) = tsb(ji, jj,jk,n1:n2-1) - tabres_child(ji,jj,jk,n1:n2-1) 126 ENDDO 127 ENDDO 128 ENDDO 129 #else 130 do jk=k1,k2 131 do jj=j1,j2 132 do ji=i1,i2 133 trbdiff(ji,jj,jk,n1:n2-1) = trb(ji,jj,jk,n1:n2) - tmask(ji,jj,jk)*tabres(ji,jj,jk,n1:n2-1)/e3t_n(ji,jj,jk) 134 enddo 135 enddo 136 enddo 137 #endif 79 138 DO jn = 1, jptra 80 139 DO jk = 1, jpkm1 -
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_top_update.F90
r6140 r8893 75 75 !! 76 76 INTEGER :: ji, jj, jk, jn 77 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child 78 REAL(wp) :: h_in(k1:k2) 79 REAL(wp) :: h_out(1:jpk) 80 INTEGER :: N_in, N_out 81 REAL(wp) :: h_diff 82 REAL(wp) :: zrho_xy 83 REAL(wp) :: tabin(k1:k2,n1:n2) 77 84 !!---------------------------------------------------------------------- 78 85 ! 79 86 IF( before ) THEN 80 ptab(i1:i2,j1:j2,k1:k2,n1:n2) = trn(i1:i2,j1:j2,k1:k2,n1:n2) 87 # if defined key_vertical 88 AGRIF_SpecialValue = -999._wp 89 zrho_xy = Agrif_rhox() * Agrif_rhoy() 90 DO jn = n1,n2-1 91 DO jk=k1,k2 92 DO jj=j1,j2 93 DO ji=i1,i2 94 ptab(ji,jj,jk,jn) = (trn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) & 95 * tmask(ji,jj,jk) * zrho_xy + (tmask(ji,jj,jk)-1)*999._wp 96 END DO 97 END DO 98 END DO 99 END DO 100 DO jk=k1,k2 101 DO jj=j1,j2 102 DO ji=i1,i2 103 ptab(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) * zrho_xy & 104 + (tmask(ji,jj,jk)-1)*999._wp 105 END DO 106 END DO 107 END DO 108 #else 109 DO jn = n1,n2-1 110 DO jk=k1,k2 111 DO jj=j1,j2 112 DO ji=i1,i2 113 tabres(ji,jj,jk,jn) = trn(ji,jj,jk,jn) 114 END DO 115 END DO 116 END DO 117 END DO 81 118 ELSE 119 #endif 120 ptab_child(:,:,:,:) = 0. 121 # if defined key_vertical 122 AGRIF_SpecialValue = 0._wp 123 DO jj=j1,j2 124 DO ji=i1,i2 125 N_in = 0 126 DO jk=k1,k2 !k2 = jpk of child grid 127 IF (ptab(ji,jj,jk,n2) == 0 ) EXIT 128 N_in = N_in + 1 129 tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)/ptab(ji,jj,jk,n2) 130 h_in(N_in) = ptab(ji,jj,jk,n2)/e1e2t(ji,jj) 131 ENDDO 132 N_out = 0 133 DO jk=1,jpk ! jpk of parent grid 134 IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF 135 N_out = N_out + 1 136 h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above 137 ENDDO 138 IF (N_in > 0) THEN !Remove this? 139 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 140 IF (h_diff < -1.e-4) THEN 141 print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out)) 142 print *,h_in(1:N_in) 143 print *,h_out(1:N_out) 144 STOP 145 ENDIF 146 DO jn=n1,n2-1 147 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) 148 ENDDO 149 ENDIF 150 ENDDO 151 ENDDO 152 #else 153 ptab_child(:,:,:,:) = ptab(:,:,:,1:jpts) 154 #endif 82 155 IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN 83 156 ! Add asselin part … … 86 159 DO jj = j1, j2 87 160 DO ji = i1, i2 88 IF( ptab (ji,jj,jk,jn) /= 0._wp ) THEN161 IF( ptab_child(ji,jj,jk,jn) /= 0._wp ) THEN 89 162 trb(ji,jj,jk,jn) = trb(ji,jj,jk,jn) & 90 & + atfp * ( ptab (ji,jj,jk,jn) &163 & + atfp * ( ptab_child(ji,jj,jk,jn) & 91 164 & - trn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 92 165 ENDIF … … 100 173 DO jj = j1, j2 101 174 DO ji = i1, i2 102 IF( ptab (ji,jj,jk,jn) /= 0._wp ) THEN103 trn(ji,jj,jk,jn) = ptab (ji,jj,jk,jn) * tmask(ji,jj,jk)175 IF( ptab_child(ji,jj,jk,jn) /= 0._wp ) THEN 176 trn(ji,jj,jk,jn) = ptab_child(ji,jj,jk,jn) * tmask(ji,jj,jk) 104 177 END IF 105 178 END DO -
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_user.F90
r8861 r8893 783 783 ! 1. Declaration of the type of variable which have to be interpolated 784 784 !--------------------------------------------------------------------- 785 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra /),trn_id)786 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra /),trn_sponge_id)785 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra+1/),trn_id) 786 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra+1/),trn_sponge_id) 787 787 788 788 ! 2. Type of interpolation
Note: See TracChangeset
for help on using the changeset viewer.
