From a8269ccfe18a58099d867c8b7df95ec983ecb788 Mon Sep 17 00:00:00 2001
From: Evgeny Mortikov <evgeny.mortikov@gmail.com>
Date: Sat, 16 Dec 2023 03:13:49 +0300
Subject: [PATCH] major code update to remove labels and goto statements
---
compile2.sh | 8 +-
drag3.f90 | 425 ++++++++++++++++++++++++--------------------------
main_drag.f90 | 8 +-
param.f90 | 29 ++--
4 files changed, 230 insertions(+), 240 deletions(-)
diff --git a/compile2.sh b/compile2.sh
index fac0309..e874355 100755
--- a/compile2.sh
+++ b/compile2.sh
@@ -1,9 +1,9 @@
#!/bin/bash
rm drag_ddt.exe *.o
-gfortran -c inputdata.f90
-gfortran -c param.f90
-gfortran -c drag3.f90
-gfortran -c main_drag.f90
+gfortran -c -Wuninitialized inputdata.f90
+gfortran -c -Wuninitialized param.f90
+gfortran -c -Wuninitialized drag3.f90
+gfortran -c -Wuninitialized main_drag.f90
gfortran -o drag_ddt.exe main_drag.o drag3.o inputdata.o param.o
diff --git a/drag3.f90 b/drag3.f90
index 41b5bd3..6b20aeb 100644
--- a/drag3.f90
+++ b/drag3.f90
@@ -1,31 +1,30 @@
- MODULE drag3
-
- USE param
- USE inputdata
+ MODULE drag3
+ USE param
+ USE inputdata
- implicit none
+ implicit none
- type, public:: data_in
- real, public:: ws, dt, st, dq, cflh, z0in
- end type
+ type, public:: data_in
+ real, public:: ws, dt, st, dq, cflh, z0in
+ end type
- type, public:: data_outdef
- real, public:: zl, ri, re, lnzuzt, zu, ztout, rith, cm, ch, ct, ckt
- end type
+ type, public:: data_outdef
+ real, public:: zl, ri, re, lnzuzt, zu, ztout, rith, cm, ch, ct, ckt
+ end type
- type, public:: data_par
- integer, public :: it=10
- end type
+ type, public:: data_par
+ integer, public :: it=10
+ end type
- type, public:: data_lutyp
- integer, public :: lu_indx=1
- end type
+ type, public:: data_lutyp
+ integer, public :: lu_indx=1
+ end type
- contains
+ contains
! SUBROUTINE surf_fluxMAS(mas1_w, mas1_dt, mas1_st, mas1_dq, out, par1, lu1)
- SUBROUTINE surf_fluxMAS(mas1_w, mas1_dt, mas1_st, mas1_dq, mas1_cflh, mas1_z0in, &
+ SUBROUTINE surf_fluxMAS(mas1_w, mas1_dt, mas1_st, mas1_dq, mas1_cflh, mas1_z0in, &
masout_zl, masout_ri, masout_re, masout_lnzuzt, masout_zu, masout_ztout,&
masout_rith, masout_cm, masout_ch, masout_ct, masout_ckt,&
par1, lu1,numst)
@@ -63,7 +62,9 @@
in%cflh=mas1_cflh(i)
in%z0in=mas1_z0in(i)
- CALL surf_flux(in, out, par1, lu1)
+ !if ((in%ws == in%ws).and.(in%dt == in%dt).and.(in%st == in%st).and.(in%dq == in%dq)) then
+ CALL surf_flux(in, out, par1, lu1)
+ !end if
masout_zl(i)=out%zl
masout_ri(i)=out%ri
@@ -80,209 +81,191 @@
END SUBROUTINE surf_fluxMAS
SUBROUTINE surf_flux(in, out, par, lu)
- type (data_in) , intent(in) :: in
- type (data_outdef) out
- type (data_par) par
- type (data_lutyp) lu
+ type (data_in) , intent(in) :: in
+ type (data_outdef) out
+ type (data_par) par
+ type (data_lutyp) lu
+
+ real ws, dt, dq, cflh, z0in
+ integer it
+ integer lu_indx
+ real zl, ri, re, lnzuzt, zu, ztin, ri_th, cm, ch, ct, ckt
+ real z0, d3, d0max, U10m, a1, y1, cimin, f, c1, u2, h0, u3, x7, d0, d00, zt, h00, ft0, an4, an5
+ real al, Tsurf, Rib, q4, t4, u, r1, f0, f4, am, o, dd, x1, y0, x0, y10, a2ch, x10, p1, p0, h, d1, f1
+ real d, c4, c1min, c0, c, b1, an0
+ integer i, j
- real ws, dt, st, dq, cflh, z0in
- integer it
- integer lu_indx
- real zl, ri, re, lnzuzt, zu, ztin, ri_th, cm, ch, ct, ckt
- real z0, d3, d0max, u1, a1, y1, cimin, f, c1, u2, h0, u3, x7, d0, d00, zt, h00, ft0, an4, an5
- real al, t1, r6, q4, t4, u, r1, f0, f4, am, o, dd, x1, y0, x0, z3, y10, a2ch, x10, p1, p0, h, d1, f1
- real d, c4, c1min, c0, c, b1, an0
- integer i, j, m
- ws=in%ws
- dt=in%dt
- st=in%st
- dq=in%dq
- cflh=in%cflh
- z0in=in%z0in
- it=par%it
+ ws=in%ws
+ dt=in%dt
+ Tsurf=in%st
+ dq=in%dq
+ cflh=in%cflh
+ z0in=in%z0in
+ it=par%it
- u=ws
- t4=dt
- t1=st
- q4=dq
- h=cflh
- z0=z0in
+ u=ws
+ t4=dt
+ q4=dq
+ h=cflh
+ z0=z0in
+
+ AN5=(Pr_t_inf_inv / Pr_t_0_inv)**4
+ D1=(2.0E0*alpha_h-AN5*alpha_m-SQRT((AN5*alpha_m)**2+4.0E0*AN5*alpha_h*(alpha_h-alpha_m)))/(2.0E0*alpha_h**2)
+ Y10=(1.0E0-alpha_m*D1)**.25E0
+ X10=(1.0E0-alpha_h*D1)**.5E0
+ P1=2.0E0*ATAN(Y10)+ALOG((Y10-1.0E0)/(Y10+1.0E0))
+ P0=ALOG((X10-1.0E0)/(X10+1.0E0))
+
+ d3=0.0e0
+ d0max=2.0e0
+
+ if(z0 < 0.0e0) then
+ !> @brief Test - definition z0 of sea surface
+ !> @details if lu_indx=2.or.lu_indx=3 call z0sea module (Ramil_Dasha)
+ d0max = 8.0e0
+ U10m=u
+ a1=0.0e0
+ y1=25.0e0
+ c1min=alog(h1/1.0e0)/kappa
+ do i=1,it
+ f=a2-2.0e0*alog(U10m)
+ do j=1,it
+ c1=(f+2.0e0*alog(y1))/kappa
+ ! looks like the check should use U10 instead of U
+ ! but note that a1 should be set = 0 in cycle beforehand
+ if(u <= 8.0e0) a1=alog(1.0e0+a3*((y1/U10m)**3))/kappa
+ c1=c1-a1
+ c1=max(c1,c1min)
+ y1=c1
+ end do
+ z0=h1*exp(-c1*kappa)
+ z0=max(z0,0.000015e0)
+ U10m=u*alog(h1/z0)/(alog(h/z0))
+ end do
+ h0=h/z0
+ u3=U10m/c1
+ else
+ ! ......parameters from viscosity sublayer......
+ h0=h/z0
+ u3=u*kappa/alog(h0)
+ end if
+
+ x7=u3*z0/an
+ if(x7 <= x8) then
+ d0=an1*alog(al1*x7)+an2
+ else
+ d0=al2*(x7**0.45e0)
+ end if
+ ! ......humidity stratification and ri-number......
+ al=g/Tsurf
+ d0=min(d0,d0max)
+ Rib=al*h*(t4+0.61e0*Tsurf*q4)/u**2
+ d00=d0
+ zt=z0/exp(d00)
+ h00=h/zt
+ ft0=alog(h00)
+ ! ......definition of r-prim......
+ an4=d1/h0
+ an5=d1/h00
+ if (abs(d0) < 1.0e-10) an5=an4
+ an5=sqrt(1.0e0-g0*an5)
+ an4=(1.0e0-alpha_m*an4)**0.25e0
+ f0=alog((x10-1.0e0)*(an5+1.0e0)/((x10+1.0e0)*(an5-1.0e0)))/Pr_t_0_inv
+ f4=2.0e0*(atan(y10)-atan(an4))+alog((y10-1.0e0)*(an4+1.0e0)/((y10+1.0e0)*(an4-1.0e0)))
+ r1=d1*f0/(f4*f4)
+ ! ......definition of dz,ta,fu,fo,fiu,fio......
+
+ if (Rib > 0.0e0) then
+ ! ......stable stratification......
+ !write (*,*) 'stable'
+
+ Rib=min(Rib,r0)
+ f=alog(h0)
+ f1=d0/f
+ a1=b4*Rib
+ a2ch=(f1+1.0e0)/Pr_t_0_inv-2.0e0*a1
+ d3=f*(sqrt(a2ch**2+4.0e0*a1*(1.0e0-a1))-a2ch)/(2.0e0*b4*(1.0e0-a1))
+ f1=b4*d3
+ f4=f+f1
+ f0=(f+d0)/Pr_t_0_inv+f1
+ o=1.0e0/Pr_t_0_inv+f1
+ am=1.0e0+f1
+ else if (Rib < r1) then
+ ! ......strong instability.....
+ !write (*,*) 'instability'
+
+ d3=d1
+ do i=1,it+1
+ d=d3/h0
+ dd=d3/h00
+ if (abs(d0) < 1.0e-10) dd=d
+ a1=(d1/d3)**(1.0e0/3.0e0)
+ x0=sqrt(1.0e0-g0*dd)
+ y0=(1.0e0-alpha_m*d)**0.25e0
+ c=alog((x0+1.0e0)/(x0-1.0e0))
+ b1=-2.0e0*atan(y0)+alog((y0+1.0e0)/(y0-1.0e0))
+ f=3.0e0*(1.0e0-a1)
+ f4=f/y10+p1+b1
+ f0=(f/x10+p0+c)/Pr_t_0_inv
+ if (i == it+1) exit
+ d3=Rib*f4**2/f0
+ end do
+ am=a1/y10
+ o=a1/(Pr_t_0_inv*x10)
+ else if (Rib > -0.001e0) then
+
+ ! ......nearly neutral......
+ write (*,*) 'neutral'
+ f4=alog(h0)
+ f0=ft0/Pr_t_0_inv
+ if (abs(d0) < 1.0e-10) f0=f4/Pr_t_0_inv
+ am=1.0e0
+ o=1.0e0/Pr_t_0_inv
+ else
+ ! ......week and semistrong instability......
+ !write (*,*) 'semistrong'
+ f1=alog(h0)
+ if (abs(d0) < 1.0e-10) ft0=f1
+ d3=Rib*Pr_t_0_inv*f1**2/ft0
+ do i=1,it+1
+ d=d3/h0
+ dd=d3/h00
+ if (abs(d0) < 1.0e-10) dd=d
+ y1=(1.0e0-alpha_m*d3)**0.25e0
+ x1=sqrt(1.0e0-g0*d3)
+ y0=(1.0e0-alpha_m*d)**0.25e0
+ x0=sqrt(1.0e0-g0*dd)
+ y0=max(y0,1.000001e0)
+ x0=max(x0,1.000001e0)
+ f4=alog((y1-1.0e0)*(y0+1.0e0)/((y1+1.0e0)*(y0-1.0e0)))+2.0e0*(atan(y1)-atan(y0))
+ f0=alog((x1-1.0e0)*(x0+1.0e0)/((x1+1.0e0)*(x0-1.0e0)))/Pr_t_0_inv
+ if (i == it + 1) exit
+ d3=Rib*f4**2/f0
+ end do
+ am=(1.0e0-alpha_m*d3)**(-0.25e0)
+ o=1.0e0/(Pr_t_0_inv*sqrt(1.0e0-g0*d3))
+ end if
+
+ ! ......computation of cu,co,k(h),alft
+ c4=kappa/f4
+ c0=kappa/f0
+ an4=kappa*c4*u*h/am
+ an0=am/o
- AN5=(A6/A0)**4
- D1=(2.0E0*G10-AN5*G4-SQRT((AN5*G4)**2+4.0E0*AN5*G10*(G10-G4)))/(2.0E0*G10**2)
- Y10=(1.0E0-G4*D1)**.25E0
- X10=(1.0E0-G10*D1)**.5E0
- P1=2.0E0*ATAN(Y10)+ALOG((Y10-1.0E0)/(Y10+1.0E0))
- P0=ALOG((X10-1.0E0)/(X10+1.0E0))
-
- d3=0.0e0
- d0max=2.0e0
+ ! ......exit......
+ out%zl=d3
+ out%ri=Rib
+ out%re=x7
+ out%lnzuzt=d00
+ out%zu=z0
+ out%ztout=zt
+ out%rith=r1
+ out%cm=c4
+ out%ch=c0
+ out%ct=an4
+ out%ckt=an0
+ return
- if(z0.lt.0.0e0) d0max=8.0e0
- if(z0.lt.0.0e0) then
- !> @brief Test - definition z0 of sea surface
- !> @details if lu_indx=2.or.lu_indx=3 call z0sea module (Ramil_Dasha)
+ END SUBROUTINE surf_flux
- u1=u
- a1=0.0e0
- y1=25.0e0
- c1min=alog(h1/1.0e0)/ap0
- do 630 i=1,it
- f=a2-2.0e0*alog(u1)
- do 570 j=1,it
- c1=(f+2.0e0*alog(y1))/ap0
- if(u.le.8.0e0) a1=alog(1.0e0+a3*((y1/u1)**3))/ap0
- c1=c1-a1
- c1=amax1(c1,c1min)
- y1=c1
- 570 continue
- z0=h1*exp(-c1*ap0)
- z0=amax1(z0,0.000015e0)
- u2=u*alog(h1/z0)/(alog(h/z0))
- u1=u2
- 630 continue
- j=1
- h0=h/z0
- u3=u1/c1
- else
-! ......parameters from viscosity sublayer......
- j=0
- h0=h/z0
- u3=u*ap0/alog(h0)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- end if
- x7=u3*z0/an
- if(x7.le.x8) then
- d0=an1*alog(al1*x7)+an2
- else
- d0=al2*(x7**0.45e0)
- end if
-! ......humidity stratification and ri-number......
- st=in%st
- al=g/t1
- d0=amin1(d0,d0max)
- r6=al*h*(t4+0.61e0*t1*q4)/u**2
- d00=d0
- zt=z0/exp(d00)
- h00=h/zt
- ft0=alog(h00)
-! ......definition of r-prim......
- an4=d1/h0
- an5=d1/h00
-!c if(d0.eq.0.0e0) an5=an4
- if (abs(d0).lt.1.0e-10) an5=an4
- an5=sqrt(1.0e0-g0*an5)
- an4=(1.0e0-g4*an4)**0.25e0
- f0=alog((x10-1.0e0)*(an5+1.0e0)/((x10+1.0e0)*(an5-1.0e0)))/a0
- f4=2.0e0*(atan(y10)-atan(an4))+alog((y10-1.0e0)*(an4+1.0e0)/((y10+1.0e0)*(an4-1.0e0)))
- r1=d1*f0/(f4*f4)
-! ......definition of dz,ta,fu,fo,fiu,fio......
- if(r6.gt.0.0e0) go to 1460
- if(r6.lt.r1) go to 1305
- if(r6.gt.-0.001e0) then
-! ......nearly neutral......
- write (*,*) 'neutral'
- f4=alog(h0)
- f0=ft0/a0
-! if(d0.eq.0.0e0) f0=f4/a0
- if (abs(d0).lt.1.0e-10) f0=f4/a0
- am=1.0e0
- o=1.0e0/a0
- go to 1570
- else
-! ......week and semistrong instability......
- write (*,*) 'semistrong'
- f1=alog(h0)
-! if(d0.eq.0.0e0)ft0=f1
- if (abs(d0).lt.1.0e-10) ft0=f1
- d3=r6*a0*f1**2/ft0
- m=1
- 1245 do 1300 i=1,it
- d=d3/h0
- dd=d3/h00
-! if(d0.eq.0.0e0)dd=d
- if (abs(d0).lt.1.0e-10) dd=d
- y1=(1.0e0-g4*d3)**0.25e0
- x1=sqrt(1.0e0-g0*d3)
- y0=(1.0e0-g4*d)**0.25e0
- x0=sqrt(1.0e0-g0*dd)
- y0=amax1(y0,1.000001e0)
- x0=amax1(x0,1.000001e0)
- f4=alog((y1-1.0e0)*(y0+1.0e0)/((y1+1.0e0)*(y0-1.0e0)))+2.0e0*(atan(y1)-atan(y0))
- f0=alog((x1-1.0e0)*(x0+1.0e0)/((x1+1.0e0)*(x0-1.0e0)))/a0
- if(m.ne.1) go to 1350
- z3=r6*f4**2/f0
- d3=z3
- 1300 continue
- m=2
- go to 1245
- 1350 am=(1.0e0-g4*d3)**(-0.25e0)
- o=1.0e0/(a0*sqrt(1.0e0-g0*d3))
- go to 1570
- end if
-! ......strong instability.....
- 1305 continue
- write (*,*) 'instability'
- d3=d1
- m=1
- 1355 do 1410 i=1,it
- d=d3/h0
- dd=d3/h00
-! if(d0.eq.0.0e0)dd=d
- if (abs(d0).lt.1.0e-10) dd=d
- a1=(d1/d3)**(1.0e0/3.0e0)
- x0=sqrt(1.0e0-g0*dd)
- y0=(1.0e0-g4*d)**0.25e0
- c=alog((x0+1.0e0)/(x0-1.0e0))
- b1=-2.0e0*atan(y0)+alog((y0+1.0e0)/(y0-1.0e0))
- f=3.0e0*(1.0e0-a1)
- f4=f/y10+p1+b1
- f0=(f/x10+p0+c)/a0
- if(m.ne.1) go to 1430
- z3=r6*f4**2/f0
- d3=z3
- 1410 continue
- m=2
- go to 1355
- 1430 am=a1/y10
- o=a1/(a0*x10)
- go to 1570
-! ......stable stratification......
- 1460 continue
- write (*,*) 'stable'
- r6=amin1(r6,r0)
- f=alog(h0)
- f1=d0/f
- a1=b4*r6
- a2ch=(f1+1.0e0)/a0-2.0e0*a1
- d3=f*(sqrt(a2ch**2+4.0e0*a1*(1.0e0-a1))-a2ch)/(2.0e0*b4*(1.0e0-a1))
- f1=b4*d3
- f4=f+f1
- f0=(f+d0)/a0+f1
- o=1.0e0/a0+f1
- am=1.0e0+f1
- 1570 continue
-! ......computation of cu,co,k(h),alft
- c4=ap0/f4
- c0=ap0/f0
- an4=ap0*c4*u*h/am
- an0=am/o
-! ......exit......
- 140 continue
- out%zl=d3
- out%ri=r6
- out%re=x7
- out%lnzuzt=d00
- out%zu=z0
- out%ztout=zt
- out%rith=r1
- out%cm=c4
- out%ch=c0
- out%ct=an4
- out%ckt=an0
-
-
-
- return
- END SUBROUTINE surf_flux
- END MODULE drag3
\ No newline at end of file
+ END MODULE drag3
\ No newline at end of file
diff --git a/main_drag.f90 b/main_drag.f90
index 6e4faff..0d18fcf 100644
--- a/main_drag.f90
+++ b/main_drag.f90
@@ -73,13 +73,13 @@
write(*,*) 'running code'
if (TEST==1) then
- filename_in='MOSAiC.txt'
+ filename_in='data/MOSAiC.txt'
filename_out='out_MOSAiC.txt'
- filename_in2='MOSAiC_zh.txt'
+ filename_in2='data/MOSAiC_zh.txt'
elseif (TEST==2) then
- filename_in='Irgason1.txt'
+ filename_in='data/Irgason1.txt'
filename_out='out_IRGASON1.txt'
- filename_in2='IRGASON_zh.txt'
+ filename_in2='data/IRGASON_zh.txt'
endif
open (1, file= filename_in, status ='old')
diff --git a/param.f90 b/param.f90
index 467f4de..16cc630 100644
--- a/param.f90
+++ b/param.f90
@@ -1,24 +1,31 @@
module param
- implicit real (a-h,o-z)
- real, parameter :: aka=.40e0
- real, parameter :: ap0=.40e0
- real, parameter :: g=9.81e0
- real, parameter :: a0=1.15e0
- real, parameter :: a6=3.5e0
- real, parameter :: g4=16.0e0
- real, parameter :: g10=16.0e0
+ implicit none
+ ! acceleration due to gravity [m/s^2]
+ real, parameter :: g = 9.81
+ ! molecular Prandtl number (air)
+ real, parameter :: Pr_m = 0.71
+ ! von Karman constant [n/d]
+ real, parameter :: kappa = 0.40
+ ! inverse Prandtl (turbulent) number in neutral conditions
+ real, parameter :: Pr_t_0_inv = 1.15
+ ! inverse Prandtl (turbulent) number in free convection
+ real, parameter :: Pr_t_inf_inv = 3.5
+
+ ! stability function coeff. [= g4 & g10 in deprecated code]
+ real, parameter :: alpha_m = 16.0
+ real, parameter :: alpha_h = 16.0
+
real, parameter :: b4=4.7e0
real, parameter :: alfam=.0144e0
real, parameter :: betam=.111e0
real, parameter :: an=.000015e0
- real, parameter :: p4=.71e0
real, parameter :: h1=10.0e0
real, parameter :: x8=16.3e0
real, parameter :: an1=5.0e0/6.0e0
real, parameter :: an2=.45e0
- real, parameter :: al1=aka*p4
+ real, parameter :: al1=kappa*Pr_m
real, parameter :: g0=1.2
- real, parameter :: al2=(.14e0*(30.0e0**an2))*(p4**.8e0)
+ real, parameter :: al2=(.14e0*(30.0e0**an2))*(Pr_m**.8e0)
real, parameter :: a2=alog(h1*(g/alfam))
real, parameter :: a3=betam*an*a2
real, parameter :: r0=.9e0/b4
--
GitLab