!! SUBROUTINE ana_initial (ng, tile, model) ! !! git $Id$ !! svn $Id: ana_initial.h 1189 2023-08-15 21:26:58Z arango $ !!====================================================================== !! Copyright (c) 2002-2023 The ROMS/TOMS Group ! !! Licensed under a MIT/X style license ! !! See License_ROMS.md ! !======================================================================= ! ! ! This subroutine sets initial conditions for momentum and tracer ! ! type variables using analytical expressions. ! ! ! !======================================================================= ! USE mod_param USE mod_grid USE mod_ncparam USE mod_ocean USE mod_stepping ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile, model ! ! Local variable declarations. ! character (len=*), parameter :: MyFile = & & __FILE__ ! #include "tile.h" ! IF (model.eq.iNLM) THEN CALL ana_NLMinitial_tile (ng, tile, model, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & GRID(ng) % h, & #ifdef SPHERICAL & GRID(ng) % lonr, & & GRID(ng) % latr, & #else & GRID(ng) % xr, & & GRID(ng) % yr, & #endif #ifdef SOLVE3D & GRID(ng) % z_r, & & OCEAN(ng) % u, & & OCEAN(ng) % v, & & OCEAN(ng) % t, & #endif & OCEAN(ng) % ubar, & & OCEAN(ng) % vbar, & & OCEAN(ng) % zeta) #ifdef TANGENT ELSE IF ((model.eq.iTLM).or.(model.eq.iRPM)) THEN CALL ana_TLMinitial_tile (ng, tile, model, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp(ng), & # ifdef SOLVE3D & nstp(ng), & & OCEAN(ng) % tl_u, & & OCEAN(ng) % tl_v, & & OCEAN(ng) % tl_t, & # endif & OCEAN(ng) % tl_ubar, & & OCEAN(ng) % tl_vbar, & & OCEAN(ng) % tl_zeta) #endif #ifdef ADJOINT ELSE IF (model.eq.iADM) THEN CALL ana_ADMinitial_tile (ng, tile, model, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & knew(ng), & # ifdef SOLVE3D & nstp(ng), & & OCEAN(ng) % ad_u, & & OCEAN(ng) % ad_v, & & OCEAN(ng) % ad_t, & # endif & OCEAN(ng) % ad_ubar, & & OCEAN(ng) % ad_vbar, & & OCEAN(ng) % ad_zeta) #endif END IF ! ! Set analytical header file name used. ! #ifdef DISTRIBUTE IF (Lanafile) THEN #else IF (Lanafile.and.(tile.eq.0)) THEN #endif ANANAME(10)=MyFile END IF ! RETURN END SUBROUTINE ana_initial ! !*********************************************************************** SUBROUTINE ana_NLMinitial_tile (ng, tile, model, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & h, & #ifdef SPHERICAL & lonr, latr, & #else & xr, yr, & #endif #ifdef SOLVE3D & z_r, & & u, v, t, & #endif & ubar, vbar, zeta) !*********************************************************************** ! USE mod_param USE mod_parallel USE mod_grid USE mod_ncparam USE mod_iounits USE mod_scalars ! #ifdef CHANNEL # ifdef DISTRIBUTE USE distribute_mod, ONLY : mp_bcasti # endif USE erf_mod, ONLY : ERF #endif USE stats_mod, ONLY : stats_2dfld #ifdef SOLVE3D USE stats_mod, ONLY : stats_3dfld #endif ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile, model integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS ! #ifdef ASSUMED_SHAPE real(r8), intent(in) :: h(LBi:,LBj:) # ifdef SPHERICAL real(r8), intent(in) :: lonr(LBi:,LBj:) real(r8), intent(in) :: latr(LBi:,LBj:) # else real(r8), intent(in) :: xr(LBi:,LBj:) real(r8), intent(in) :: yr(LBi:,LBj:) # endif # ifdef SOLVE3D real(r8), intent(in) :: z_r(LBi:,LBj:,:) real(r8), intent(out) :: u(LBi:,LBj:,:,:) real(r8), intent(out) :: v(LBi:,LBj:,:,:) real(r8), intent(out) :: t(LBi:,LBj:,:,:,:) # endif real(r8), intent(out) :: ubar(LBi:,LBj:,:) real(r8), intent(out) :: vbar(LBi:,LBj:,:) real(r8), intent(out) :: zeta(LBi:,LBj:,:) #else # ifdef SPHERICAL real(r8), intent(in) :: lonr(LBi:UBi,LBj:UBj) real(r8), intent(in) :: latr(LBi:UBi,LBj:UBj) # else real(r8), intent(in) :: xr(LBi:UBi,LBj:UBj) real(r8), intent(in) :: yr(LBi:UBi,LBj:UBj) # endif real(r8), intent(in) :: h(LBi:UBi,LBj:UBj) # ifdef SOLVE3D real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng)) real(r8), intent(out) :: u(LBi:UBi,LBj:UBj,N(ng),2) real(r8), intent(out) :: v(LBi:UBi,LBj:UBj,N(ng),2) real(r8), intent(out) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng)) # endif real(r8), intent(out) :: ubar(LBi:UBi,LBj:UBj,:) real(r8), intent(out) :: vbar(LBi:UBi,LBj:UBj,:) real(r8), intent(out) :: zeta(LBi:UBi,LBj:UBj,:) #endif ! ! Local variable declarations. ! logical, save :: first = .TRUE. ! integer :: Iless, Iplus, i, itrc, j, k ! #ifdef CHANNEL real(r8), parameter :: guscale = 40.0E+03_r8 real(r8), parameter :: u0 = 1.6_r8 #endif real(r8) :: depth, dx, val1, val2, val3, val4, x, x0, y, y0 ! TYPE (T_STATS), save :: Stats(7) ! ubar, vbar, zeta, u, v, t, s #include "set_bounds.h" ! !----------------------------------------------------------------------- ! Initialize field statistics structure. !----------------------------------------------------------------------- ! IF (first) THEN first=.FALSE. DO i=1,SIZE(Stats,1) Stats(i) % checksum=0_i8b Stats(i) % count=0 Stats(i) % min=Large Stats(i) % max=-Large Stats(i) % avg=0.0_r8 Stats(i) % rms=0.0_r8 END DO END IF ! !----------------------------------------------------------------------- ! Initial conditions for 2D momentum (m/s) components. !----------------------------------------------------------------------- ! #if defined CHANNEL && !defined ONLY_TS_IC y0=0.5_r8*el(ng) DO j=JstrT,JendT DO i=IstrP,IendT val1=(GRID(ng)%yu(i,j)-y0)/guscale ubar(i,j,1)=u0*EXP(-val1*val1)/6.0_r8 END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT vbar(i,j,1)=0.0_r8 END DO END DO #elif defined SOLITON x0=2.0_r8*xl(ng)/3.0_r8 y0=0.5_r8*el(ng) val1=0.395_r8 val2=0.771_r8*(val1*val1) IF (ng.eq.1) THEN DO j=JstrT,JendT DO i=IstrP,IendT x=0.5_r8*(xr(i-1,j)+xr(i,j))-x0 y=0.5_r8*(yr(i-1,j)+yr(i,j))-y0 val3=EXP(-val1*x) val4=val2*((2.0_r8*val3/(1.0_r8+(val3*val3)))**2) ubar(i,j,1)=0.25_r8*val4*(6.0_r8*y*y-9.0_r8)* & & EXP(-0.5_r8*y*y) END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT x=0.5_r8*(xr(i,j-1)+xr(i,j))-x0 y=0.5_r8*(yr(i,j-1)+yr(i,j))-y0 val3=EXP(-val1*x) val4=val2*((2.0_r8*val3/(1.0_r8+(val3*val3)))**2) vbar(i,j,1)=2.0_r8*val4*y*(-2.0_r8*val1*TANH(val1*x))* & & EXP(-0.5_r8*y*y) END DO END DO ELSE DO j=JstrT,JendT DO i=IstrP,IendT ubar(i,j,1)=0.0_r8 END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT vbar(i,j,1)=0.0_r8 END DO END DO END IF #elif defined RIVERPLUME2 DO j=JstrT,JendT DO i=IstrP,IendT ubar(i,j,1)=0.0_r8 END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT vbar(i,j,1)=-0.05_r8 END DO END DO #elif defined SED_TEST1 val1=100.0_r8 DO j=JstrT,JendT DO i=IstrP,IendT ubar(i,j,1)=-10.0_r8/(10.0_r8+9.0E-06_r8*REAL(i,r8)*val1) END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT vbar(i,j,1)=0.0_r8 END DO END DO #elif defined SED_TOY DO j=JstrT,JendT DO i=IstrP,IendT ubar(i,j,1)=0.0_r8 END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT vbar(i,j,1)=0.0_r8 END DO END DO #elif defined TEST_CHAN val1=100.0_r8 DO j=JstrT,JendT DO i=IstrP,IendT ubar(i,j,1)=0.0_r8 END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT vbar(i,j,1)=0.0_r8 END DO END DO #else DO j=JstrT,JendT DO i=IstrP,IendT ubar(i,j,1)=0.0_r8 END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT vbar(i,j,1)=0.0_r8 END DO END DO #endif ! ! Report statistics. ! CALL stats_2dfld (ng, tile, iNLM, u2dvar, Stats(1), & & LBi, UBi, LBj, UBj, ubar(:,:,1)) IF (DOMAIN(ng)%NorthEast_Corner(tile)) THEN WRITE (stdout,10) TRIM(Vname(2,idUbar))//': '// & & TRIM(Vname(1,idUbar)), & & ng, Stats(1)%min, Stats(1)%max END IF CALL stats_2dfld (ng, tile, iNLM, v2dvar, Stats(2), & & LBi, UBi, LBj, UBj, vbar(:,:,1)) IF (DOMAIN(ng)%NorthEast_Corner(tile)) THEN WRITE (stdout,10) TRIM(Vname(2,idVbar))//': '// & & TRIM(Vname(1,idVbar)), & & ng, Stats(2)%min, Stats(2)%max END IF ! !----------------------------------------------------------------------- ! Initial conditions for free-surface (m). !----------------------------------------------------------------------- ! #if defined CHANNEL && !defined ONLY_TS_IC y0=0.5_r8*el(ng) # ifdef SOLVE3D DO j=JstrT,JendT DO i=IstrT,IendT val1=(yr(i,j)-y0)/guscale val2=-u0*guscale*GRID(ng)%f(i,j)*SQRT(pi)/(12.0_r8*g) zeta(i,j,1)=val2*ERF(val1) END DO END DO # else DO j=JstrT,JendT DO i=IstrT,IendT val1=(yr(i,j)-y0)/guscale val2=-0.5_r8*u0*guscale*GRID(ng)%f(i,j)*sqrt(pi)/g zeta(i,j,1)=val2*ERF(val1) END DO END DO # endif # ifdef DISTRIBUTE CALL mp_bcasti (ng, model, exit_flag) ! in case of error in ERF # endif #elif defined KELVIN !! val1=1.0_r8 ! zeta0 !! val2=2.0_r8*pi/(12.42_r8*3600.0_r8) ! M2 Tide period DO j=JstrT,JendT DO i=IstrT,IendT !! zeta(i,j,1)=val1* & !! & EXP(-GRID(ng)%f(i,j)*GRID(ng)%yp(i,j)/ & !! & SQRT(g*GRID(ng)%h(i,j)))* & !! & COS(val2*GRID(ng)%xp(i,j)/ & !! & SQRT(g*GRID(ng)%h(i,j))) zeta(i,j,1)=0.0_r8 END DO END DO #elif defined SOLITON IF (ng.eq.1) THEN x0=2.0_r8*xl(ng)/3.0_r8 y0=0.5_r8*el(ng) val1=0.395_r8 val2=0.771_r8*(val1*val1) DO j=JstrT,JendT DO i=IstrT,IendT x=xr(i,j)-x0 y=yr(i,j)-y0 val3=EXP(-val1*x) val4=val2*((2.0_r8*val3/(1.0_r8+(val3*val3)))**2) zeta(i,j,1)=0.25_r8*val4*(6.0_r8*y*y+3.0_r8)* & & EXP(-0.5_r8*y*y) END DO END DO ELSE DO j=JstrT,JendT DO i=IstrT,IendT zeta(i,j,1)=0.0_r8 END DO END DO END IF #elif defined SED_TEST1 val1=100.0_r8 DO j=JstrT,JendT DO i=IstrT,IendT zeta(i,j,1)=9.0E-06_r8*REAL(i,r8)*val1 END DO END DO #elif defined TEST_CHAN DO j=JstrT,JendT DO i=IstrT,IendT zeta(i,j,1)=-0.4040_r8*REAL(i,r8)/REAL(Lm(ng)+1,r8) END DO END DO #else DO j=JstrT,JendT DO i=IstrT,IendT zeta(i,j,1)=0.0_r8 END DO END DO #endif ! ! Report statistics. ! CALL stats_2dfld (ng, tile, iNLM, r2dvar, Stats(3), & & LBi, UBi, LBj, UBj, zeta(:,:,1)) IF (DOMAIN(ng)%NorthEast_Corner(tile)) THEN WRITE (stdout,10) TRIM(Vname(2,idFsur))//': '// & & TRIM(Vname(1,idFsur)), & & ng, Stats(3)%min, Stats(3)%max END IF #ifdef SOLVE3D ! !----------------------------------------------------------------------- ! Initial conditions for 3D momentum components (m/s). !----------------------------------------------------------------------- ! # if defined CHANNEL && !defined ONLY_TS_IC y0=0.5_r8*el(ng) DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrP,IendT val1=(GRID(ng)%yu(i,j)-y0)/guscale val2=(z_r(i,j,k)+z_r(i-1,j,k))/(h(i,j)+h(i-1,j)) val3=u0*(0.5_r8+val2+(0.5_r8*val2*val2))*EXP(-val1*val1) u(i,j,k,1)=val3 END DO END DO END DO DO k=1,N(ng) DO j=JstrP,JendT DO i=IstrT,IendT v(i,j,k,1)=0.0_r8 END DO END DO END DO # elif defined RIVERPLUME2 DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrP,IendT u(i,j,k,1)=0.0_r8 END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT v(i,j,k,1)=-0.05_r8*LOG((h(i,j)+z_r(i,j,k))/Zob(ng))/ & & (LOG(h(i,j)/Zob(ng))-1.0_r8+Zob(ng)/h(i,j)) END DO END DO END DO # elif defined SED_TEST1 DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrP,IendT u(i,j,k,1)=-1.0_r8*LOG((h(i,j)+z_r(i,j,k))/Zob(ng))/ & & (LOG(h(i,j)/Zob(ng))-1.0_r8+Zob(ng)/h(i,j)) END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT v(i,j,k,1)=0.0_r8 END DO END DO END DO # elif defined SED_TOY DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrP,IendT u(i,j,k,1)=1.0_r8 !! u(i,j,k,1)=-1.0_r8 !! u(i,j,k,1)=0.0_r8 END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT v(i,j,k,1)=0.0_r8 END DO END DO END DO # else DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrP,IendT u(i,j,k,1)=0.0_r8 END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT v(i,j,k,1)=0.0_r8 END DO END DO END DO # endif ! ! Report statistics. ! CALL stats_3dfld (ng, tile, iNLM, u3dvar, Stats(4), & & LBi, UBi, LBj, UBj, 1, N(ng), u(:,:,:,1)) IF (DOMAIN(ng)%NorthEast_Corner(tile)) THEN WRITE (stdout,10) TRIM(Vname(2,idUvel))//': '// & & TRIM(Vname(1,idUvel)), & & ng, Stats(4)%min, Stats(4)%max END IF CALL stats_3dfld (ng, tile, iNLM, v3dvar, Stats(5), & & LBi, UBi, LBj, UBj, 1, N(ng), v(:,:,:,1)) IF (DOMAIN(ng)%NorthEast_Corner(tile)) THEN WRITE (stdout,10) TRIM(Vname(2,idVvel))//': '// & & TRIM(Vname(1,idVvel)), & & ng, Stats(5)%min, Stats(5)%max END IF ! !----------------------------------------------------------------------- ! Initial conditions for tracer type variables. !----------------------------------------------------------------------- ! ! Set initial conditions for potential temperature (Celsius) and ! salinity (PSU). ! # if defined BENCHMARK val1=(44.69_r8/39.382_r8)**2 val2=val1*(rho0*800.0_r8/g)*(5.0E-05_r8/((42.689_r8/44.69_r8)**2)) DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=val2*EXP(z_r(i,j,k)/800.0_r8)* & & (0.6_r8-0.4_r8*TANH(z_r(i,j,k)/800.0_r8)) # ifdef SALINITY t(i,j,k,1,isalt)=35.0_r8 # endif END DO END DO END DO # elif defined BASIN val1=(44.69_r8/39.382_r8)**2 val2=val1*(rho0*800.0_r8/g)*(5.0E-05_r8/((42.689_r8/44.69_r8)**2)) DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=val2*EXP(z_r(i,j,k)/800.0_r8)* & & (0.6_r8-0.4_r8*TANH(z_r(i,j,k)/800.0_r8)) END DO END DO END DO # elif defined BL_TEST DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT val1=TANH(1.1_r8*z_r(i,j,k)+11.0_r8) t(i,j,k,1,itemp)=T0(ng)+6.25_r8*val1 # ifdef SALINITY t(i,j,k,1,isalt)=S0(ng)-0.75_r8*val1 # endif END DO END DO END DO # elif defined CHANNEL y0=0.5_r8*el(ng) DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT val1=(yr(i,j)-y0)/guscale val2=-0.5_r8*u0*guscale*GRID(ng)%f(i,j)*SQRT(pi)/ & & (Tcoef(ng)*g*h(i,j)) val3=(val2*ERF(val1)+T0(ng))*(1.0_r8+z_r(i,j,k)/h(i,j)) t(i,j,k,1,itemp)=val3 END DO END DO END DO # ifdef DISTRIBUTE CALL mp_bcasti (ng, model, exit_flag) ! in case of error in ERF # endif # elif defined CANYON DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=3.488_r8*EXP(z_r(i,j,k)/800.0_r8)* & & (1.0_r8-(2.0_r8/3.0_r8)* & & TANH(z_r(i,j,k)/800.0_r8)) END DO END DO END DO # elif defined CHANNEL_NECK DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=20.0_r8 !! t(i,j,k,1,itemp)=14.0_r8+8.0_r8*EXP(z_r(i,j,k)/50.0_r8) END DO END DO END DO # elif defined COUPLING_TEST val1=40.0_r8 DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=val1*EXP(z_r(i,j,k)/800.0_r8)* & & (0.6_r8-0.4_r8*TANH(z_r(i,j,k)/800.0_r8))+ & & 1.5_r8 # ifdef SALINITY t(i,j,k,1,isalt)=35.0_r8 # endif END DO END DO END DO # elif defined DOUBLE_GYRE val1=(44.69_r8/39.382_r8)**2 val2=val1*(rho0*100.0_r8/g)*(5.0E-5_r8/((42.689_r8/44.69_r8)**2)) DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT val3=T0(ng)+val2*EXP(z_r(i,j,k)/100.0_r8)* & & (10.0_r8-0.4_r8*TANH(z_r(i,j,k)/100.0_r8)) val4=yr(i,j)/el(ng) t(i,j,k,1,itemp)=val3-3.0_r8*val4 # ifdef SALINITY t(i,j,k,1,isalt)=34.5_r8-0.001_r8*z_r(i,j,k)-val4 # endif END DO END DO END DO # elif defined ESTUARY_TEST DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=10.0_r8 # ifdef SALINITY IF (xr(i,j).le.30000.0_r8) then t(i,j,k,1,isalt)=30.0_r8 ELSEIF (xr(i,j).le.80000.0_r8) then t(i,j,k,1,isalt)=(80000.0_r8-xr(i,j))/50000.0_r8*30.0_r8 ELSE t(i,j,k,1,isalt)=0.0_r8 END IF # endif END DO END DO END DO # elif defined FLT_TEST DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=T0(ng) END DO END DO END DO # elif defined GRAV_ADJ DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,MIN((Lm(ng)+1)/2,IendT) t(i,j,k,1,itemp)=T0(ng)+5.0_r8 # ifdef SALINITY t(i,j,k,1,isalt)=0.0_r8 # endif END DO DO i=MAX((Lm(ng)+1)/2+1,IstrT),IendT t(i,j,k,1,itemp)=T0(ng) # ifdef SALINITY t(i,j,k,1,isalt)=S0(ng) # endif END DO !! DO i=IstrT,IendT !! IF (i.lt.Lm(ng)/2) THEN !! t(i,j,k,1,itemp)=T0(ng)+5.0_r8 !! ELSE IF (i.eq.Lm(ng)/2) THEN !! t(i,j,k,1,itemp)=T0(ng)+4.0_r8 !! ELSE IF (i.eq.Lm(ng)/2+1) THEN !! t(i,j,k,1,itemp)=T0(ng)+1.0_r8 !! ELSE !! t(i,j,k,1,itemp)=T0(ng) !! END IF !! END DO END DO END DO # elif defined LAB_CANYON DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=-659.34183_r8*z_r(i,j,k) END DO END DO END DO # elif defined LAKE_SIGNELL DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=10.0_r8 # ifdef SALINITY t(i,j,k,1,isalt)=30.0_r8 # endif END DO END DO END DO # elif defined LMD_TEST DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=MIN(13.0_r8, & & 7.0_r8+0.2_r8*(z_r(i,j,k)+50.0_r8)) # ifdef SALINITY t(i,j,k,1,isalt)=35.0_r8 # endif END DO END DO END DO # elif defined MIXED_LAYER DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=10.0_r8+3.0_r8*(z_r(i,j,k)+h(i,j))/ & & h(i,j) # ifdef SALINITY t(i,j,k,1,isalt)=S0(ng) # endif END DO END DO END DO # elif defined NJ_BIGHT DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT depth=z_r(i,j,k) IF (depth.ge.-15.0_r8) THEN t(i,j,k,1,itemp)= 2.049264257728403E+01_r8-depth* & & (2.640850848793918E-01_r8+depth* & & (2.751125328535212E-01_r8+depth* & & (9.207489761648872E-02_r8+depth* & & (1.449075725742839E-02_r8+depth* & & (1.078215685912076E-03_r8+depth* & & (3.240318053903974E-05_r8+ & & 1.262826857690271E-07_r8*depth)))))) # ifdef SALINITY t(i,j,k,1,isalt)= 3.066489149193135E+01_r8-depth* & & (1.476725262946735E-01_r8+depth* & & (1.126455760313399E-01_r8+depth* & & (3.900923281871022E-02_r8+depth* & & (6.939014937447098E-03_r8+depth* & & (6.604436696792939E-04_r8+depth* & & (3.191792361954220E-05_r8+ & & 6.177352634409320E-07_r8*depth)))))) # endif ELSE t(i,j,k,1,itemp)=14.6_r8+ & & 6.70_r8*TANH(1.1_r8*depth+15.9_r8) # ifdef SALINITY t(i,j,k,1,isalt)=31.3_r8- & & 0.55_r8*TANH(1.1_r8*depth+15.9_r8) # endif END IF END DO END DO END DO # elif defined OVERFLOW DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=T0(ng)-0.5_r8*T0(ng)*(1.0_r8+ & & TANH((yr(i,j)-60000.0_r8)/2000.0_r8)) END DO END DO END DO # elif defined RIVERPLUME1 DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=T0(ng)+0.01_r8*REAL(k,r8) # ifdef SALINITY t(i,j,k,1,isalt)=S0(ng) # endif END DO END DO END DO # elif defined RIVERPLUME2 DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=T0(ng) # ifdef SALINITY t(i,j,k,1,isalt)=S0(ng) # endif END DO END DO END DO # elif defined SEAMOUNT DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=T0(ng)+7.5_r8*EXP(z_r(i,j,k)/1000.0_r8) END DO END DO END DO # elif defined SED_TEST1 DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=20.0_r8 # ifdef SALINITY t(i,j,k,1,isalt)=0.0_r8 # endif END DO END DO END DO # elif defined UPWELLING DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=T0(ng)+8.0_r8*EXP(z_r(i,j,k)/50.0_r8) !! t(i,j,k,1,itemp)=T0(ng)+(z_r(i,j,k)+75.0_r8)/150.0_r8+ !! & 4.0_r8*(1.0_r8+TANH((z_r(i,j,k)+35.0_r8)/ !! & 6.5_r8)) # ifdef SALINITY !! t(i,j,k,1,isalt)=1.0E-04_r8*yr(i,j)-S0(ng) t(i,j,k,1,isalt)=S0(ng) !! IF (j.lt.Mm(ng)/2) THEN !! t(i,j,k,1,isalt)=0.0_r8 !! ELSE IF (j.eq.Mm(ng)/2) THEN !! t(i,j,k,1,isalt)=0.5_r8 !! ELSE IF (j.gt.Mm(ng)/2) THEN !! t(i,j,k,1,isalt)=1.0_r8 !! END IF # endif END DO END DO END DO # elif defined WINDBASIN DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=20.0_r8 ! homogeneous !! t(i,j,k,1,itemp)=14.0_r8+8.0_r8*EXP(z_r(i,j,k)/50.0_r8)- & !! & T0(ng) ! stratified END DO END DO END DO # else DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT t(i,j,k,1,itemp)=T0(ng) # ifdef SALINITY t(i,j,k,1,isalt)=S0(ng) # endif END DO END DO END DO # endif ! ! Report statistics. ! DO itrc=1,NAT CALL stats_3dfld (ng, tile, iNLM, u3dvar, Stats(itrc+5), & & LBi, UBi, LBj, UBj, 1, N(ng), t(:,:,:,1,itrc)) IF (DOMAIN(ng)%NorthEast_Corner(tile)) THEN WRITE (stdout,10) TRIM(Vname(2,idTvar(itrc)))//': '// & & TRIM(Vname(1,idTvar(itrc))), & & ng, Stats(itrc+5)%min, Stats(itrc+5)%max END IF END DO #endif ! 10 FORMAT (3x,' ANA_INITIAL - ',a,/,19x, & & '(Grid = ',i2.2,', Min = ',1p,e15.8,0p, & & ' Max = ',1p,e15.8,0p,')') ! RETURN END SUBROUTINE ana_NLMinitial_tile #ifdef TANGENT ! !*********************************************************************** SUBROUTINE ana_TLMinitial_tile (ng, tile, model, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & kstp, & # ifdef SOLVE3D & nstp, & & tl_u, tl_v, tl_t, & # endif & tl_ubar, tl_vbar, tl_zeta) !*********************************************************************** ! USE mod_param USE mod_scalars ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile, model integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: kstp # ifdef SOLVE3D integer, intent(in) :: nstp # endif ! # ifdef ASSUMED_SHAPE # ifdef SOLVE3D real(r8), intent(out) :: tl_u(LBi:,LBj:,:,:) real(r8), intent(out) :: tl_v(LBi:,LBj:,:,:) real(r8), intent(out) :: tl_t(LBi:,LBj:,:,:,:) # endif real(r8), intent(out) :: tl_ubar(LBi:,LBj:,:) real(r8), intent(out) :: tl_vbar(LBi:,LBj:,:) real(r8), intent(out) :: tl_zeta(LBi:,LBj:,:) # else # ifdef SOLVE3D real(r8), intent(out) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2) real(r8), intent(out) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2) real(r8), intent(out) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng)) # endif real(r8), intent(out) :: tl_ubar(LBi:UBi,LBj:UBj,:) real(r8), intent(out) :: tl_vbar(LBi:UBi,LBj:UBj,:) real(r8), intent(out) :: tl_zeta(LBi:UBi,LBj:UBj,:) # endif ! ! Local variable declarations. ! integer :: i, itrc, j, k # include "set_bounds.h" ! !----------------------------------------------------------------------- ! Initial conditions for tangent linear 2D momentum (s/m) components. !----------------------------------------------------------------------- ! DO j=JstrT,JendT DO i=IstrP,IendT tl_ubar(i,j,kstp)=0.0_r8 END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT tl_vbar(i,j,kstp)=0.0_r8 END DO END DO ! !----------------------------------------------------------------------- ! Initial conditions for tangent linear free-surface (1/m). !----------------------------------------------------------------------- ! DO j=JstrT,JendT DO i=IstrT,IendT tl_zeta(i,j,kstp)=0.0_r8 END DO END DO # ifdef SOLVE3D ! !----------------------------------------------------------------------- ! Initial conditions for tangent linear 3D momentum components (s/m). !----------------------------------------------------------------------- ! DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrP,IendT tl_u(i,j,k,nstp)=0.0_r8 END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT tl_v(i,j,k,nstp)=0.0_r8 END DO END DO END DO ! !----------------------------------------------------------------------- ! Initial conditions for tangent linear active tracers (1/Tunits). !----------------------------------------------------------------------- ! DO itrc=1,NAT DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT tl_t(i,j,k,nstp,itrc)=0.0_r8 END DO END DO END DO END DO # endif ! RETURN END SUBROUTINE ana_TLMinitial_tile #endif #ifdef ADJOINT ! !*********************************************************************** SUBROUTINE ana_ADMinitial_tile (ng, tile, model, & & LBi, UBi, LBj, UBj, & & IminS, ImaxS, JminS, JmaxS, & & knew, & # ifdef SOLVE3D & nstp, & & ad_u, ad_v, ad_t, & # endif & ad_ubar, ad_vbar, ad_zeta) !*********************************************************************** ! USE mod_param USE mod_scalars ! ! Imported variable declarations. ! integer, intent(in) :: ng, tile, model integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: IminS, ImaxS, JminS, JmaxS integer, intent(in) :: knew # ifdef SOLVE3D integer, intent(in) :: nstp # endif ! # ifdef ASSUMED_SHAPE # ifdef SOLVE3D real(r8), intent(out) :: ad_u(LBi:,LBj:,:,:) real(r8), intent(out) :: ad_v(LBi:,LBj:,:,:) real(r8), intent(out) :: ad_t(LBi:,LBj:,:,:,:) # endif real(r8), intent(out) :: ad_ubar(LBi:,LBj:,:) real(r8), intent(out) :: ad_vbar(LBi:,LBj:,:) real(r8), intent(out) :: ad_zeta(LBi:,LBj:,:) # else # ifdef SOLVE3D real(r8), intent(out) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2) real(r8), intent(out) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2) real(r8), intent(out) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng)) # endif real(r8), intent(out) :: ad_ubar(LBi:UBi,LBj:UBj,:) real(r8), intent(out) :: ad_vbar(LBi:UBi,LBj:UBj,:) real(r8), intent(out) :: ad_zeta(LBi:UBi,LBj:UBj,:) # endif ! ! Local variable declarations. ! integer :: i, itrc, j, k # include "set_bounds.h" ! !----------------------------------------------------------------------- ! Initial conditions for adjoint 2D momentum (s/m) components. !----------------------------------------------------------------------- ! DO j=JstrT,JendT DO i=IstrP,IendT ad_ubar(i,j,knew)=0.0_r8 END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT ad_vbar(i,j,knew)=0.0_r8 END DO END DO ! !----------------------------------------------------------------------- ! Initial conditions for adjoint free-surface (1/m). !----------------------------------------------------------------------- ! DO j=JstrT,JendT DO i=IstrT,IendT ad_zeta(i,j,knew)=0.0_r8 END DO END DO # ifdef SOLVE3D ! !----------------------------------------------------------------------- ! Initial conditions for adjoint 3D momentum components (s/m). !----------------------------------------------------------------------- ! DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrP,IendT ad_u(i,j,k,nstp)=0.0_r8 END DO END DO DO j=JstrP,JendT DO i=IstrT,IendT ad_v(i,j,k,nstp)=0.0_r8 END DO END DO END DO ! !----------------------------------------------------------------------- ! Initial conditions for adjoint active tracers (1/Tunits). !----------------------------------------------------------------------- ! DO itrc=1,NAT DO k=1,N(ng) DO j=JstrT,JendT DO i=IstrT,IendT ad_t(i,j,k,nstp,itrc)=0.0_r8 END DO END DO END DO END DO # endif ! RETURN END SUBROUTINE ana_ADMinitial_tile #endif