[782] | 1 | MODULE agrif_oce |
---|
[1605] | 2 | !!====================================================================== |
---|
[782] | 3 | !! *** MODULE agrif_oce *** |
---|
[1605] | 4 | !! AGRIF : define in memory AGRIF variables |
---|
[782] | 5 | !!---------------------------------------------------------------------- |
---|
[1605] | 6 | !! History : 2.0 ! 2007-12 (R. Benshila) Original code |
---|
[782] | 7 | !!---------------------------------------------------------------------- |
---|
[1605] | 8 | #if defined key_agrif |
---|
[782] | 9 | !!---------------------------------------------------------------------- |
---|
[1605] | 10 | !! 'key_agrif' AGRIF zoom |
---|
| 11 | !!---------------------------------------------------------------------- |
---|
[782] | 12 | USE par_oce ! ocean parameters |
---|
| 13 | USE dom_oce ! domain parameters |
---|
[5656] | 14 | |
---|
[782] | 15 | IMPLICIT NONE |
---|
[2715] | 16 | PRIVATE |
---|
[782] | 17 | |
---|
[2715] | 18 | PUBLIC agrif_oce_alloc ! routine called by nemo_init in nemogcm.F90 |
---|
[9031] | 19 | #if defined key_vertical |
---|
| 20 | PUBLIC reconstructandremap ! remapping routine |
---|
| 21 | #endif |
---|
[1605] | 22 | ! !!* Namelist namagrif: AGRIF parameters |
---|
[11542] | 23 | LOGICAL , PUBLIC :: ln_agrif_2way = .TRUE. !: activate two way nesting |
---|
| 24 | LOGICAL , PUBLIC :: ln_spc_dyn = .FALSE. !: use zeros (.false.) or not (.true.) in |
---|
| 25 | !: bdys dynamical fields interpolation |
---|
[5656] | 26 | REAL(wp), PUBLIC :: rn_sponge_tra = 2800. !: sponge coeff. for tracers |
---|
| 27 | REAL(wp), PUBLIC :: rn_sponge_dyn = 2800. !: sponge coeff. for dynamics |
---|
| 28 | LOGICAL , PUBLIC :: ln_chk_bathy = .FALSE. !: check of parent bathymetry |
---|
[11244] | 29 | ! |
---|
| 30 | INTEGER , PUBLIC, PARAMETER :: nn_sponge_len = 2 !: Sponge width (in number of parent grid points) |
---|
[5656] | 31 | LOGICAL , PUBLIC :: spongedoneT = .FALSE. !: tracer sponge layer indicator |
---|
| 32 | LOGICAL , PUBLIC :: spongedoneU = .FALSE. !: dynamics sponge layer indicator |
---|
| 33 | LOGICAL , PUBLIC :: lk_agrif_fstep = .TRUE. !: if true: first step |
---|
| 34 | LOGICAL , PUBLIC :: lk_agrif_debug = .FALSE. !: if true: print debugging info |
---|
[1605] | 35 | |
---|
[5656] | 36 | LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_tsn |
---|
| 37 | # if defined key_top |
---|
| 38 | LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_trn |
---|
| 39 | # endif |
---|
| 40 | LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_u |
---|
| 41 | LOGICAL , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tabspongedone_v |
---|
[11219] | 42 | INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: utint_stage |
---|
| 43 | INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vtint_stage |
---|
[9019] | 44 | REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fsaht_spu, fsaht_spv !: sponge diffusivities |
---|
| 45 | REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fsahm_spt, fsahm_spf !: sponge viscosities |
---|
[2715] | 46 | |
---|
[9019] | 47 | ! Barotropic arrays used to store open boundary data during time-splitting loop: |
---|
[11205] | 48 | REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ubdy, vbdy, hbdy |
---|
[5656] | 49 | |
---|
[9019] | 50 | |
---|
| 51 | INTEGER, PUBLIC :: tsn_id ! AGRIF profile for tracers interpolation and update |
---|
| 52 | INTEGER, PUBLIC :: un_interp_id, vn_interp_id ! AGRIF profiles for interpolations |
---|
| 53 | INTEGER, PUBLIC :: un_update_id, vn_update_id ! AGRIF profiles for udpates |
---|
| 54 | INTEGER, PUBLIC :: tsn_sponge_id, un_sponge_id, vn_sponge_id ! AGRIF profiles for sponge layers |
---|
[5656] | 55 | # if defined key_top |
---|
[9019] | 56 | INTEGER, PUBLIC :: trn_id, trn_sponge_id |
---|
[5656] | 57 | # endif |
---|
[9019] | 58 | INTEGER, PUBLIC :: unb_id, vnb_id, ub2b_interp_id, vb2b_interp_id |
---|
| 59 | INTEGER, PUBLIC :: ub2b_update_id, vb2b_update_id |
---|
| 60 | INTEGER, PUBLIC :: e3t_id, e1u_id, e2v_id, sshn_id |
---|
| 61 | INTEGER, PUBLIC :: scales_t_id |
---|
| 62 | INTEGER, PUBLIC :: avt_id, avm_id, en_id ! TKE related identificators |
---|
| 63 | INTEGER, PUBLIC :: umsk_id, vmsk_id |
---|
| 64 | INTEGER, PUBLIC :: kindic_agr |
---|
| 65 | |
---|
[1605] | 66 | !!---------------------------------------------------------------------- |
---|
[9598] | 67 | !! NEMO/NST 4.0 , NEMO Consortium (2018) |
---|
[1605] | 68 | !! $Id$ |
---|
[10068] | 69 | !! Software governed by the CeCILL license (see ./LICENSE) |
---|
[2715] | 70 | !!---------------------------------------------------------------------- |
---|
| 71 | CONTAINS |
---|
| 72 | |
---|
| 73 | INTEGER FUNCTION agrif_oce_alloc() |
---|
| 74 | !!---------------------------------------------------------------------- |
---|
| 75 | !! *** FUNCTION agrif_oce_alloc *** |
---|
| 76 | !!---------------------------------------------------------------------- |
---|
[5656] | 77 | INTEGER, DIMENSION(2) :: ierr |
---|
| 78 | !!---------------------------------------------------------------------- |
---|
| 79 | ierr(:) = 0 |
---|
| 80 | ! |
---|
[11219] | 81 | ALLOCATE( fsaht_spu(jpi,jpj), fsaht_spv(jpi,jpj), & |
---|
| 82 | & fsahm_spt(jpi,jpj), fsahm_spf(jpi,jpj), & |
---|
| 83 | & tabspongedone_tsn(jpi,jpj), & |
---|
| 84 | & utint_stage(jpi,jpj), vtint_stage(jpi,jpj), & |
---|
[5656] | 85 | # if defined key_top |
---|
| 86 | & tabspongedone_trn(jpi,jpj), & |
---|
| 87 | # endif |
---|
| 88 | & tabspongedone_u (jpi,jpj), & |
---|
| 89 | & tabspongedone_v (jpi,jpj), STAT = ierr(1) ) |
---|
| 90 | |
---|
[11205] | 91 | ALLOCATE( ubdy(jpi,jpj), vbdy(jpi,jpj), hbdy(jpi,jpj), STAT = ierr(2) ) |
---|
[5656] | 92 | |
---|
| 93 | agrif_oce_alloc = MAXVAL(ierr) |
---|
| 94 | ! |
---|
[2715] | 95 | END FUNCTION agrif_oce_alloc |
---|
| 96 | |
---|
[9031] | 97 | #if defined key_vertical |
---|
| 98 | SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout) |
---|
| 99 | !!---------------------------------------------------------------------- |
---|
| 100 | !! *** FUNCTION reconstructandremap *** |
---|
| 101 | !!---------------------------------------------------------------------- |
---|
| 102 | IMPLICIT NONE |
---|
| 103 | INTEGER N, Nout |
---|
| 104 | REAL(wp) tabin(N), tabout(Nout) |
---|
| 105 | REAL(wp) hin(N), hout(Nout) |
---|
| 106 | REAL(wp) coeffremap(N,3),zwork(N,3) |
---|
| 107 | REAL(wp) zwork2(N+1,3) |
---|
| 108 | INTEGER jk |
---|
| 109 | DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8 |
---|
| 110 | REAL(wp) q,q01,q02,q001,q002,q0 |
---|
| 111 | REAL(wp) z_win(1:N+1), z_wout(1:Nout+1) |
---|
| 112 | REAL(wp),PARAMETER :: dpthin = 1.D-3 |
---|
| 113 | INTEGER :: k1, kbox, ktop, ka, kbot |
---|
| 114 | REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop |
---|
| 115 | |
---|
| 116 | z_win(1)=0.; z_wout(1)= 0. |
---|
| 117 | DO jk=1,N |
---|
| 118 | z_win(jk+1)=z_win(jk)+hin(jk) |
---|
| 119 | ENDDO |
---|
| 120 | |
---|
| 121 | DO jk=1,Nout |
---|
| 122 | z_wout(jk+1)=z_wout(jk)+hout(jk) |
---|
| 123 | ENDDO |
---|
| 124 | |
---|
| 125 | DO jk=2,N |
---|
| 126 | zwork(jk,1)=1./(hin(jk-1)+hin(jk)) |
---|
| 127 | ENDDO |
---|
| 128 | |
---|
| 129 | DO jk=2,N-1 |
---|
| 130 | q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1)) |
---|
| 131 | zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1) |
---|
| 132 | zwork(jk,3)=q0 |
---|
| 133 | ENDDO |
---|
| 134 | |
---|
| 135 | DO jk= 2,N |
---|
| 136 | zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1)) |
---|
| 137 | ENDDO |
---|
| 138 | |
---|
| 139 | coeffremap(:,1) = tabin(:) |
---|
| 140 | |
---|
| 141 | DO jk=2,N-1 |
---|
| 142 | q001 = hin(jk)*zwork2(jk+1,1) |
---|
| 143 | q002 = hin(jk)*zwork2(jk,1) |
---|
| 144 | IF (q001*q002 < 0) then |
---|
| 145 | q001 = 0. |
---|
| 146 | q002 = 0. |
---|
| 147 | ENDIF |
---|
| 148 | q=zwork(jk,2) |
---|
| 149 | q01=q*zwork2(jk+1,1) |
---|
| 150 | q02=q*zwork2(jk,1) |
---|
| 151 | IF (abs(q001) > abs(q02)) q001 = q02 |
---|
| 152 | IF (abs(q002) > abs(q01)) q002 = q01 |
---|
| 153 | |
---|
| 154 | q=(q001-q002)*zwork(jk,3) |
---|
| 155 | q001=q001-q*hin(jk+1) |
---|
| 156 | q002=q002+q*hin(jk-1) |
---|
| 157 | |
---|
| 158 | coeffremap(jk,3)=coeffremap(jk,1)+q001 |
---|
| 159 | coeffremap(jk,2)=coeffremap(jk,1)-q002 |
---|
| 160 | |
---|
| 161 | zwork2(jk,1)=(2.*q001-q002)**2 |
---|
| 162 | zwork2(jk,2)=(2.*q002-q001)**2 |
---|
| 163 | ENDDO |
---|
| 164 | |
---|
| 165 | DO jk=1,N |
---|
| 166 | IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN |
---|
| 167 | coeffremap(jk,3) = coeffremap(jk,1) |
---|
| 168 | coeffremap(jk,2) = coeffremap(jk,1) |
---|
| 169 | zwork2(jk,1) = 0. |
---|
| 170 | zwork2(jk,2) = 0. |
---|
| 171 | ENDIF |
---|
| 172 | ENDDO |
---|
| 173 | |
---|
| 174 | DO jk=2,N |
---|
| 175 | q002=max(zwork2(jk-1,2),dsmll) |
---|
| 176 | q001=max(zwork2(jk,1),dsmll) |
---|
| 177 | zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002) |
---|
| 178 | ENDDO |
---|
| 179 | |
---|
| 180 | zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3) |
---|
| 181 | zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3) |
---|
| 182 | |
---|
| 183 | DO jk=1,N |
---|
| 184 | q01=zwork2(jk+1,3)-coeffremap(jk,1) |
---|
| 185 | q02=coeffremap(jk,1)-zwork2(jk,3) |
---|
| 186 | q001=2.*q01 |
---|
| 187 | q002=2.*q02 |
---|
| 188 | IF (q01*q02<0) then |
---|
| 189 | q01=0. |
---|
| 190 | q02=0. |
---|
| 191 | ELSEIF (abs(q01)>abs(q002)) then |
---|
| 192 | q01=q002 |
---|
| 193 | ELSEIF (abs(q02)>abs(q001)) then |
---|
| 194 | q02=q001 |
---|
| 195 | ENDIF |
---|
| 196 | coeffremap(jk,2)=coeffremap(jk,1)-q02 |
---|
| 197 | coeffremap(jk,3)=coeffremap(jk,1)+q01 |
---|
| 198 | ENDDO |
---|
| 199 | |
---|
| 200 | zbot=0.0 |
---|
| 201 | kbot=1 |
---|
| 202 | DO jk=1,Nout |
---|
| 203 | ztop=zbot !top is bottom of previous layer |
---|
| 204 | ktop=kbot |
---|
| 205 | IF (ztop.GE.z_win(ktop+1)) then |
---|
| 206 | ktop=ktop+1 |
---|
| 207 | ENDIF |
---|
| 208 | |
---|
| 209 | zbot=z_wout(jk+1) |
---|
| 210 | zthk=zbot-ztop |
---|
| 211 | |
---|
| 212 | IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN |
---|
| 213 | |
---|
| 214 | kbot=ktop |
---|
| 215 | DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N) |
---|
| 216 | kbot=kbot+1 |
---|
| 217 | ENDDO |
---|
| 218 | zbox=zbot |
---|
| 219 | DO k1= jk+1,Nout |
---|
| 220 | IF (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN |
---|
| 221 | exit !thick layer |
---|
| 222 | ELSE |
---|
| 223 | zbox=z_wout(k1+1) !include thin adjacent layers |
---|
| 224 | IF(zbox.EQ.z_wout(Nout+1)) THEN |
---|
| 225 | exit !at bottom |
---|
| 226 | ENDIF |
---|
| 227 | ENDIF |
---|
| 228 | ENDDO |
---|
| 229 | zthk=zbox-ztop |
---|
| 230 | |
---|
| 231 | kbox=ktop |
---|
| 232 | DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N) |
---|
| 233 | kbox=kbox+1 |
---|
| 234 | ENDDO |
---|
| 235 | |
---|
| 236 | IF(ktop.EQ.kbox) THEN |
---|
| 237 | IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN |
---|
| 238 | IF(hin(kbox).GT.dpthin) THEN |
---|
| 239 | q001 = (zbox-z_win(kbox))/hin(kbox) |
---|
| 240 | q002 = (ztop-z_win(kbox))/hin(kbox) |
---|
| 241 | q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002) |
---|
| 242 | q02=q01-1.+(q001+q002) |
---|
| 243 | q0=1.-q01-q02 |
---|
| 244 | ELSE |
---|
| 245 | q0 = 1.0 |
---|
| 246 | q01 = 0. |
---|
| 247 | q02 = 0. |
---|
| 248 | ENDIF |
---|
| 249 | tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3) |
---|
| 250 | ELSE |
---|
| 251 | tabout(jk) = tabin(kbox) |
---|
| 252 | ENDIF |
---|
| 253 | ELSE |
---|
| 254 | IF(ktop.LE.jk .AND. kbox.GE.jk) THEN |
---|
| 255 | ka = jk |
---|
| 256 | ELSEIF (kbox-ktop.GE.3) THEN |
---|
| 257 | ka = (kbox+ktop)/2 |
---|
| 258 | ELSEIF (hin(ktop).GE.hin(kbox)) THEN |
---|
| 259 | ka = ktop |
---|
| 260 | ELSE |
---|
| 261 | ka = kbox |
---|
| 262 | ENDIF !choose ka |
---|
| 263 | |
---|
| 264 | offset=coeffremap(ka,1) |
---|
| 265 | |
---|
| 266 | qtop = z_win(ktop+1)-ztop !partial layer thickness |
---|
| 267 | IF(hin(ktop).GT.dpthin) THEN |
---|
| 268 | q=(ztop-z_win(ktop))/hin(ktop) |
---|
| 269 | q01=q*(q-1.) |
---|
| 270 | q02=q01+q |
---|
| 271 | q0=1-q01-q02 |
---|
| 272 | ELSE |
---|
| 273 | q0 = 1. |
---|
| 274 | q01 = 0. |
---|
| 275 | q02 = 0. |
---|
| 276 | ENDIF |
---|
| 277 | |
---|
| 278 | tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop |
---|
| 279 | |
---|
| 280 | DO k1= ktop+1,kbox-1 |
---|
| 281 | tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1) |
---|
| 282 | ENDDO !k1 |
---|
| 283 | |
---|
| 284 | qbot = zbox-z_win(kbox) !partial layer thickness |
---|
| 285 | IF(hin(kbox).GT.dpthin) THEN |
---|
| 286 | q=qbot/hin(kbox) |
---|
| 287 | q01=(q-1.)**2 |
---|
| 288 | q02=q01-1.+q |
---|
| 289 | q0=1-q01-q02 |
---|
| 290 | ELSE |
---|
| 291 | q0 = 1.0 |
---|
| 292 | q01 = 0. |
---|
| 293 | q02 = 0. |
---|
| 294 | ENDIF |
---|
| 295 | |
---|
| 296 | tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot |
---|
| 297 | |
---|
| 298 | rpsum=1.0d0/zthk |
---|
| 299 | tabout(jk)=offset+tsum*rpsum |
---|
| 300 | |
---|
| 301 | ENDIF !single or multiple layers |
---|
| 302 | ELSE |
---|
| 303 | IF (jk==1) THEN |
---|
| 304 | write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1) |
---|
| 305 | ENDIF |
---|
| 306 | tabout(jk) = tabout(jk-1) |
---|
| 307 | |
---|
| 308 | ENDIF !normal:thin layer |
---|
| 309 | ENDDO !jk |
---|
| 310 | |
---|
| 311 | return |
---|
| 312 | end subroutine reconstructandremap |
---|
[2715] | 313 | #endif |
---|
[9031] | 314 | |
---|
| 315 | #endif |
---|
[1605] | 316 | !!====================================================================== |
---|
[782] | 317 | END MODULE agrif_oce |
---|