Skip to content

Commit 5001c53

Browse files
mjuckermjucker
mjucker
authored and
mjucker
committed
corrected T_fact,tau_fact bugs
Two bugs, one concerning the variable name T_fact (and not Te_fact), and one concerning tau_strat have been corrected.
1 parent 921fdfd commit 5001c53

File tree

1 file changed

+39
-40
lines changed

1 file changed

+39
-40
lines changed

‎fms_riga_hs_jucker/src/atmos_param/hs_forcing/hs_forcing.f90

Lines changed: 39 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -74,7 +74,7 @@ module hs_forcing_mod
7474
real :: local_heating_latwidth=0.4 ! radians latitude Used only when local_heating_option='Gaussian'
7575
real :: local_heating_sigwidth=0.11 ! sigma height Used only when local_heating_option='Gaussian'
7676
real :: local_heating_sigcenter=0.3 ! sigma height Used only when local_heating_option='Gaussian'
77-
logical :: polar_heating_option=.false. ! want to add some heating over the pole?
77+
logical :: polar_heating_option=.false. ! want to add some heating over the pole?
7878
real :: polar_heating_srfamp=0.0 ! Degrees per day Used only when polar heating_option='true'
7979
real :: polar_heating_latwidth=0.0 ! radians latitude Used only when polar_heating_option='true'
8080
real :: polar_heating_latcenter=0.0 ! radians latitude Used only when polar_heating_option='true'
@@ -150,7 +150,7 @@ module hs_forcing_mod
150150
equilibrium_tau_option,equilibrium_tau_file, & !mj
151151
p_hs,p_bd,A_NH_0,A_NH_1,A_SH_0,A_SH_1,A_s, & !mj
152152
phi_N,phi_S,tau_t,tau_N_p,tau_S_p,delta_phi, & !mj
153-
Te_fact,tau_fact, & !mj
153+
T_fact,tau_fact, & !mj
154154
tau_m,do_seasonal_cycle,days_per_year !mj
155155

156156
!-----------------------------------------------------------------------
@@ -276,7 +276,7 @@ subroutine hs_forcing ( is, ie, js, je, dt, Time, lon, lat, p_half, p_full, &
276276
call get_number_tracers(MODEL_ATMOS, num_tracers=num_tracers)
277277
n_hum = get_tracer_index(MODEL_ATMOS, 'sphum') !mj
278278
n_pv = get_tracer_index(MODEL_ATMOS, 'pv') !mj
279-
279+
280280
if(num_tracers == size(rdt,4)) then
281281
do n = 1, size(rdt,4)
282282
flux = trflux
@@ -293,12 +293,12 @@ subroutine hs_forcing ( is, ie, js, je, dt, Time, lon, lat, p_half, p_full, &
293293
elseif(n == get_tracer_index(MODEL_ATMOS,'methane')) then !mj methane tracer
294294
rst = rm(:,:,:,n) + dt*rdt(:,:,:,n) !mj
295295
call methane_source_sink ( flux, sink, p_half, rst, rtnd, kbot ) !mj
296-
rdt(:,:,:,n) = rdt(:,:,:,n) + rtnd !mj
296+
rdt(:,:,:,n) = rdt(:,:,:,n) + rtnd !mj
297297
elseif(n == n_pv) then
298298
rst = rm(:,:,:,n) + dt*rdt(:,:,:,n) !mj
299299
tst = tm + dt*tdt !mj
300300
call pv_tracer(vom, tst, lat, p_full, dt, rst, rtnd)
301-
rdt(:,:,:,n) = rdt(:,:,:,n) + rtnd !mj
301+
rdt(:,:,:,n) = rdt(:,:,:,n) + rtnd !mj
302302
elseif(n == get_tracer_index(MODEL_ATMOS, 'APV')) then !mj Blocking index as Schwierz2004
303303
rst = rm(:,:,:,n) + dt*rdt(:,:,:,n) !mj
304304
tst = tm + dt*tdt !mj
@@ -309,7 +309,7 @@ subroutine hs_forcing ( is, ie, js, je, dt, Time, lon, lat, p_half, p_full, &
309309
pvt = rm(:,:,:,n_pv) + dt*rdt(:,:,:,n_pv) !mj
310310
call apv_tracer(pvt,p_full,dt,rst,rtnd)
311311
endif
312-
rdt(:,:,:,n) = rdt(:,:,:,n) + rtnd !mj
312+
rdt(:,:,:,n) = rdt(:,:,:,n) + rtnd !mj
313313
else !mj
314314
if (query_method('tracer_sms', MODEL_ATMOS, n, scheme, params)) then
315315
if (uppercase(trim(scheme)) == 'NONE') cycle
@@ -347,7 +347,7 @@ subroutine hs_forcing_init ( axes, Time, lonb, latb )
347347
integer, intent(in) :: axes(4)
348348
type(time_type), intent(in) :: Time
349349
real, intent(in), optional, dimension(:,:) :: lonb, latb
350-
350+
351351

352352
!-----------------------------------------------------------------------
353353
integer unit, io, ierr
@@ -405,7 +405,7 @@ subroutine hs_forcing_init ( axes, Time, lonb, latb )
405405

406406
! If positive, damping time units are (1/s), value is the inverse of damping time.
407407
! If negative, damping time units are (days), value is the damping time. It is converted to (1/s)
408-
408+
409409
if (ka < 0.) then
410410
tka = -1./(86400*ka)
411411
else
@@ -488,15 +488,15 @@ subroutine hs_forcing_init ( axes, Time, lonb, latb )
488488
&trim(equilibrium_tau_option) == 'strat_file' ) then
489489
call interpolator_init (tau_interp, trim(equilibrium_tau_file)//'.nc', lonb, latb, data_out_of_bounds=(/CONSTANT/),vert_interp=(/INTERP_LINEAR_P/))
490490
endif
491-
491+
492492

493493
module_is_initialized = .true.
494494

495495
end subroutine hs_forcing_init
496496

497497
!#######################################################################
498498

499-
subroutine hs_forcing_end
499+
subroutine hs_forcing_end
500500

501501
!-----------------------------------------------------------------------
502502
!
@@ -524,7 +524,7 @@ subroutine hs_forcing_end
524524
call interpolator_end(u_interp)
525525
call interpolator_end(v_interp)
526526
endif
527-
527+
528528
module_is_initialized = .false.
529529

530530
end subroutine hs_forcing_end
@@ -634,7 +634,7 @@ subroutine newtonian_damping ( Time, lat, ps, p_full, p_half, t, tdt, teq, tau,
634634

635635
t_star(:,:) = t_zero - delh*sin_lat_2(:,:) - eps_sc*sin_lat(:,:)
636636
if ( .not. pv_sat_flag) then
637-
tstr (:,:) = t_strat
637+
tstr (:,:) = t_strat
638638
else
639639
tstr (:,:) = t_tropopause
640640
endif
@@ -688,7 +688,7 @@ subroutine newtonian_damping ( Time, lat, ps, p_full, p_half, t, tdt, teq, tau,
688688
endwhere
689689
else
690690
call error_mesg ('hs_forcing_nml', &
691-
'"'//trim(equilibrium_tau_option)//'" is not a valid value for equilibrium_tau_option',FATAL)
691+
'"'//trim(equilibrium_tau_option)//'" is not a valid value for equilibrium_tau_option',FATAL)
692692
endif
693693

694694

@@ -707,7 +707,7 @@ subroutine newtonian_damping ( Time, lat, ps, p_full, p_half, t, tdt, teq, tau,
707707
p_norm(:,:) = p_full(:,:,k)/pref
708708

709709
do l = 1,sat_levs - 1
710-
where (p_norm(:,:) < p_tropopause .and. p_norm(:,:) > sat_p(l+1) .and. p_norm(:,:) <= sat_p(l))
710+
where (p_norm(:,:) < p_tropopause .and. p_norm(:,:) > sat_p(l+1) .and. p_norm(:,:) <= sat_p(l))
711711
t_sat(:,:) = sat_t(l) * (p_norm(:,:)/sat_p(l))**(-rdgas*sat_g(l)/grav);
712712
end where
713713
end do
@@ -734,7 +734,7 @@ subroutine newtonian_damping ( Time, lat, ps, p_full, p_half, t, tdt, teq, tau,
734734
phipi_S = phi_S*pif
735735
where ( lat .ge. phipi_S .and. lat .le. phipi_N )
736736
P3 = -1.96e-9*(lat/pif)**4 - 1.15e-5*(lat/pif)**2 + 1
737-
elsewhere ( lat .lt. phipi_S )
737+
elsewhere ( lat .lt. phipi_S )
738738
P3 = -1.96e-9*(phi_S)**4 - 1.15e-5*(phi_S)**2 + 1
739739
elsewhere ( lat .gt. phipi_N )
740740
P3 = -1.96e-9*(phi_N)**4 - 1.15e-5*(phi_N)**2 + 1
@@ -756,7 +756,7 @@ subroutine newtonian_damping ( Time, lat, ps, p_full, p_half, t, tdt, teq, tau,
756756
P4 = (P3 - 1.)*p_norm + 1.
757757
elsewhere ( p_full(:,:,k) <= p_1 )
758758
P4 = P3
759-
endwhere
759+
endwhere
760760
! Add polar amplitudes
761761
p_norm(:,:) = log(p_full(:,:,k)/p_t)/(log(p_1/1000.e2) - log(p_t/1000.e2))
762762
p_n = log(p_1 /p_t)/(log(p_1/1000.e2) - log(p_t/1000.e2))
@@ -811,7 +811,7 @@ subroutine newtonian_damping ( Time, lat, ps, p_full, p_half, t, tdt, teq, tau,
811811
Pp = ( logphs - log(p_full) )/( logphs - logpbd )
812812
teq = Pp*teq_strat + (1. - Pp)*teq
813813
elsewhere( p_full .lt. p_bd )
814-
teq = teq_strat
814+
teq = teq_strat
815815
end where
816816
end if
817817
!! same for damping rate
@@ -844,9 +844,9 @@ subroutine newtonian_damping ( Time, lat, ps, p_full, p_half, t, tdt, teq, tau,
844844
do k=1,size(t,3)
845845
P4 = min(1.,( Pp(:,:,k) - P2 )/P3)
846846
tau_strat(:,:,k) = ( tau_strat(:,:,k) - tau_m )*P4 + tau_m
847+
! scale profile
848+
tau_strat(:,:,k) = (1.-tau_fact)*tau_t + tau_fact*tau_strat(:,:,k)
847849
enddo
848-
! scale profile
849-
tau_strat(:,:,k) = (1.-tau_fact)*tau_t + tau_fact*tau_strat(:,:,k)
850850
! convert days to seconds
851851
tau_strat = tau_strat*86400
852852
! convert to damping rate
@@ -863,7 +863,7 @@ subroutine newtonian_damping ( Time, lat, ps, p_full, p_half, t, tdt, teq, tau,
863863
Pp = ( logphs - log(p_full) )/( logphs - logpbd )
864864
tdamp = Pp*tau_strat + (1. - Pp)*tdamp
865865
elsewhere( p_full .lt. p_bd )
866-
tdamp = tau_strat
866+
tdamp = tau_strat
867867
end where
868868
endif
869869

@@ -1006,7 +1006,7 @@ subroutine sphum_source_sink ( flux, damp, p_full, r, rdt, s_geo, t, dt, kbot )
10061006
if (rdamp > 0.) rdamp = 1./rdamp
10071007

10081008
!------------ surface source and global sink ---------------------------
1009-
1009+
10101010
source(:,:,:)=0.0
10111011
sink(:,:,:)=0.0 !mj
10121012
sea_surf=0 !mj
@@ -1029,7 +1029,7 @@ subroutine sphum_source_sink ( flux, damp, p_full, r, rdt, s_geo, t, dt, kbot )
10291029
kb = size(r,3)
10301030
qsat = RDGAS*ES0*exp(-HLV*(1./t - 1./TFREEZE)/RVGAS)/RVGAS !mj
10311031
qsat = qsat/p_full !mj
1032-
source(:,:,kb) = max(0.,-vkf*(r(:,:,kb)-qsat(:,:,kb))) !mj
1032+
source(:,:,kb) = max(0.,-vkf*(r(:,:,kb)-qsat(:,:,kb))) !mj
10331033
source(:,:,kb) = source(:,:,kb)*sea_surf(:,:) !mj mountains
10341034
endif
10351035

@@ -1039,9 +1039,9 @@ subroutine sphum_source_sink ( flux, damp, p_full, r, rdt, s_geo, t, dt, kbot )
10391039
! sink = 1. + qsat*HLV*HLV/(RVGAS*CP_AIR*t*t) !mj as in Frierson 2006 JAS
10401040
! sink = (r - qsat)/sink/dt !mj continued
10411041
end where
1042-
1042+
10431043
! sink = sink + rdamp*r !mj add global sink
1044-
1044+
10451045
! sink = 0.
10461046
! sink = r/dt
10471047
! source = 0.
@@ -1066,13 +1066,13 @@ subroutine age_source_sink ( flux, damp, p_half, r, rdt, kbot )
10661066
integer :: i, j, kb
10671067
real :: rdamp
10681068
!-----------------------------------------------------------------------
1069-
1069+
10701070
rdamp = damp
10711071
if (rdamp < 0.) rdamp = -86400.*rdamp ! convert days to seconds
10721072
if (rdamp > 0.) rdamp = 1./rdamp
10731073

10741074
!------------ simple surface source and no sink --------------------
1075-
1075+
10761076
source=0.0
10771077
sink =0.0
10781078

@@ -1092,7 +1092,7 @@ subroutine age_source_sink ( flux, damp, p_half, r, rdt, kbot )
10921092

10931093
sink = rdamp*r
10941094

1095-
rdt = source - sink
1095+
rdt = source - sink
10961096

10971097
!-----------------------------------------------------------------------
10981098

@@ -1118,7 +1118,7 @@ subroutine methane_source_sink ( flux, damp, p_half, r, rdt, kbot )
11181118
if (rdamp > 0.) rdamp = 1./rdamp
11191119

11201120
!------------ surface source and global sink ---------------------------
1121-
1121+
11221122
source(:,:,:)=0.0
11231123
sink(:,:,:)=0.0 !mj
11241124

@@ -1134,14 +1134,14 @@ subroutine methane_source_sink ( flux, damp, p_half, r, rdt, kbot )
11341134
kb = size(r,3)
11351135
source(:,:,kb) = -vkf*(r(:,:,kb)-flux) !mj tracer value fixed to trflux at bottom
11361136
endif
1137-
1137+
11381138

11391139
sink = rdamp*r
1140-
1140+
11411141

11421142
rdt = source - sink
11431143

1144-
1144+
11451145

11461146
!-----------------------------------------------------------------------
11471147

@@ -1160,8 +1160,8 @@ subroutine pv_tracer ( vorn, temp, lat, p_full, dt, r, rdt )
11601160
! real,dimension(size(rdt,1),size(rdt,2)) :: wa,dTheta
11611161
real :: wa,dTheta
11621162
!-----------------------------------------------------------------------
1163-
1164-
1163+
1164+
11651165
do k=2,size(rdt,3) !actual value
11661166
do j=1,size(rdt,2)
11671167
do i=1,size(rdt,1)
@@ -1173,7 +1173,7 @@ subroutine pv_tracer ( vorn, temp, lat, p_full, dt, r, rdt )
11731173
enddo
11741174
rdt(:,:,1)=rdt(:,:,2)
11751175
!!$ rdt = vorn
1176-
1176+
11771177
if(dt.gt.0.0)then ! normal pv tracer
11781178
rdt = (rdt-r)/dt
11791179
endif
@@ -1197,8 +1197,8 @@ subroutine apv_tracer ( pv, p_half, dt, r, rdt )
11971197
p_high = 50000.
11981198
del_p = 0.
11991199
apv = 0.
1200-
1201-
do k=2,size(rdt,3)
1200+
1201+
do k=2,size(rdt,3)
12021202
do j=1,size(rdt,2)
12031203
do i=1,size(rdt,1)
12041204
if(p_half(i,j,k) >= p_low .and. p_half(i,j,k+1) <= p_high)then
@@ -1215,7 +1215,7 @@ subroutine apv_tracer ( pv, p_half, dt, r, rdt )
12151215
enddo
12161216
rdt(:,:,1)=rdt(:,:,2)
12171217
!!$ rdt = vorn
1218-
1218+
12191219
rdt = (rdt-r)/dt
12201220

12211221
!-----------------------------------------------------------------------
@@ -1338,8 +1338,8 @@ subroutine local_heating ( Time, is, js, lon, lat, ps, p_full, p_half, surf_geop
13381338
if(z_full(i,j,k)<=12e3) then
13391339
z_factor = sin((pi*z_full(i,j,k))/12e3)
13401340
tdt(i,j,k) = srfamp*lat_factor(i,j)*z_factor
1341-
else
1342-
tdt(i,j,k) = 0
1341+
else
1342+
tdt(i,j,k) = 0
13431343
endif
13441344
enddo
13451345
enddo
@@ -1456,4 +1456,3 @@ end subroutine get_tau
14561456
!#######################################################################
14571457

14581458
end module hs_forcing_mod
1459-

0 commit comments

Comments
 (0)