diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..b2e5ffe858c28defec07efbbc2c0ba0bb753e9fc
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1 @@
+/drag.exe
diff --git a/compile2.sh b/compile2.sh
old mode 100644
new mode 100755
index e1433f88e7ddfcdc6627204ad266a985b7e85a89..fac0309cab857b44bbe7059823b359fcdaf67da3
--- a/compile2.sh
+++ b/compile2.sh
@@ -3,7 +3,6 @@
rm drag_ddt.exe *.o
gfortran -c inputdata.f90
gfortran -c param.f90
-gfortran -c prmt.f90
gfortran -c drag3.f90
gfortran -c 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 039a999c735642db86d2f47cbcd2b1711d7b5715..41b5bd398b4da3d0637263e9252bbd00ac555567 100644
--- a/drag3.f90
+++ b/drag3.f90
@@ -2,9 +2,10 @@
USE param
USE inputdata
-
- !implicit real (a-h, o-z)
+
implicit none
+
+
type, public:: data_in
real, public:: ws, dt, st, dq, cflh, z0in
end type
@@ -16,20 +17,80 @@
type, public:: data_par
integer, public :: it=10
end type
-
+
+ type, public:: data_lutyp
+ integer, public :: lu_indx=1
+ end type
contains
-
- SUBROUTINE surf_flux(in, out, par)
+
+ ! 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, &
+ 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)
+
+ real, dimension (numst) :: mas1_w
+ real, dimension (numst) :: mas1_dt
+ real, dimension (numst) :: mas1_st
+ real, dimension (numst) :: mas1_dq
+ real, dimension (numst) :: mas1_cflh
+ real, dimension (numst) :: mas1_z0in
+
+ real, dimension (numst) :: masout_zl
+ real, dimension (numst) :: masout_ri
+ real, dimension (numst) :: masout_re
+ real, dimension (numst) :: masout_lnzuzt
+ real, dimension (numst) :: masout_zu
+ real, dimension (numst) :: masout_ztout
+ real, dimension (numst) :: masout_rith
+ real, dimension (numst) :: masout_cm
+ real, dimension (numst) :: masout_ch
+ real, dimension (numst) :: masout_ct
+ real, dimension (numst) :: masout_ckt
+
+ type (data_in) in
+ type (data_outdef) out
+ type (data_par) par1
+ type (data_lutyp) lu1
+ integer i,numst
+
+ do i=1,numst
+ in%ws=mas1_w(i)
+ in%dt=mas1_dt(i)
+ in%st=mas1_st(i)
+ in%dq=mas1_dq(i)
+ in%cflh=mas1_cflh(i)
+ in%z0in=mas1_z0in(i)
+
+ CALL surf_flux(in, out, par1, lu1)
+
+ masout_zl(i)=out%zl
+ masout_ri(i)=out%ri
+ masout_re(i)=out%re
+ masout_lnzuzt(i)=out%lnzuzt
+ masout_zu(i)=out%zu
+ masout_ztout(i)=out%ztout
+ masout_rith(i)=out%rith
+ masout_cm(i)=out%cm
+ masout_ch(i)=out%ch
+ masout_ct(i)=out%ct
+ masout_ckt(i)=out%ckt
+ enddo
+ 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_par) par
+ type (data_lutyp) lu
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, h1, ap0, f, a2, c1, u2, h0, u3, x7, x8, an1, an2, d0, d00, zt, h00, ft0, an4, an5
- real al, al2, an, g, t1, r6, q4, t4, u, g0, r1, f0, f4, a0, am, o, dd, x1, y0, x0, z3, y10, a2ch, x10, p1, p0, h, d1, f1
+ 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
@@ -46,44 +107,48 @@
q4=dq
h=cflh
z0=z0in
+
+ 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
- !=data_in%ws
- !4=data_in%dt
- !4=data_in%dq
- !=data_in%cflh
- !0=data_in%z0
+
if(z0.lt.0.0e0) d0max=8.0e0
if(z0.lt.0.0e0) then
-! ......definition z0 of sea surface......
- !call z0sea module (ramil_dasha)
- !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
+ !> @brief Test - definition z0 of sea surface
+ !> @details if lu_indx=2.or.lu_indx=3 call z0sea module (Ramil_Dasha)
+
+ 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)
+ j=0
+ h0=h/z0
+ u3=u*ap0/alog(h0)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end if
x7=u3*z0/an
diff --git a/inputdata.f90 b/inputdata.f90
index 8fb673770fd57abb890deac91524726ffb6f2157..3543420ad45a569feac67f279a641116dde66c7b 100644
--- a/inputdata.f90
+++ b/inputdata.f90
@@ -1,7 +1,14 @@
module INPUTDATA
REAl, DIMENSION (6) :: AR1
REAl, DIMENSION (11) :: AR2
- INTEGER, PARAMETER :: IT =1
+ INTEGER, PARAMETER :: TEST=1 ! probably IT before renaming
+ INTEGER nums,ioer,mk
+ REAL HFX, MFX, zL, betta
+ REAL U, T4,c0,c4, T1,H,z0h
+ REAL ws, deltaT, semisumT, D00, Z0, ZT, deltaQ
+ REAL R6,R1
+ REAL AN5, D1, Y10, X10, P1, P0
+
!C*====================================================================
!C* .....DEFENITION OF DRAG AND HEAT EXCHANGE COEFFICIENTS...... =
!C* DETAILS OF ALGORITM ARE GIVEN IN: =
diff --git a/main_drag.f90 b/main_drag.f90
index 6c809aff54932813bc9aa783c8e3c3ac0d0358bf..6e4faff854ec8cef011bea8841888cd258ef01f1 100644
--- a/main_drag.f90
+++ b/main_drag.f90
@@ -7,27 +7,164 @@
type (data_in):: data_in1
type (data_outdef) :: data_outdef1
-
-
+
type (data_par) :: data_par1
-
- open (1,file='4_ddt.txt')
-
- do i = 1,1000000
- read (1,*,end=100) data_in1%ws, data_in1%dt
-
-
- CALL surf_flux(data_in1, data_outdef1, data_par1)
+ type (data_lutyp) :: data_lutyp1
+ integer :: numst, i
+ real :: cflh, z0in
+ character(len = 50) :: filename_in
+ character(len = 50) :: filename_out
+ character(len = 50) :: filename_in2
+
+ type :: datatype_inMAS1
+ real, allocatable :: mas_w(:) !
+ real, allocatable :: mas_dt(:)
+ real, allocatable :: mas_st(:)
+ real, allocatable :: mas_dq(:)
+ real, allocatable :: mas_cflh(:)
+ real, allocatable :: mas_z0in(:)
+ end type
+ type(datatype_inMAS1) :: data_inMAS
+
+ !input
+ ! mas_w - abs(wind velocity) at constant flux layer (cfl) hight (m/s)
+ ! mas_dt - difference between potential temperature at cfl hight and at surface ( deg. k)
+ ! mas_st - semi-sum of potential temperature at cfl hight and and at surface ( deg. k)
+ ! mas_dq - difference between humidity at cfl hight and a surface ( gr/gr )
+ ! cflh - - cfl hight ( m )
+ ! z0in=0.01 - roughness of surface ( m );
+ ! it - number of iterations
+ ! lu_indx - 1 for land, 2 for sea, 3 for lake
+ ! test - file input
+
+ type :: datatype_outMAS1
+ real, allocatable :: masout_zl(:)
+ real, allocatable :: masout_ri(:)
+ real, allocatable :: masout_re(:)
+ real, allocatable :: masout_lnzuzt(:)
+ real, allocatable :: masout_zu(:)
+ real, allocatable :: masout_ztout(:)
+ real, allocatable :: masout_rith(:)
+ real, allocatable :: masout_cm(:)
+ real, allocatable :: masout_ch(:)
+ real, allocatable :: masout_ct(:)
+ real, allocatable :: masout_ckt(:)
+ end type
+
+ type(datatype_outMAS1) :: data_outMAS
+
+ !output
+ !masout_zl - non-dimensional cfl hight
+ !masout_ri - richardson number
+ !masout_re - reynods number
+ !masout_lnzuzt - ln(zu/zt)
+ !masout_zu - dynamical roughness zu (m)
+ !masout_ztout - thermal roughness zt (m)
+ !masout_rith - critical richardson number
+ !masout_cm - transfer coefficient for momentum
+ !masout_ch - transfer coefficient fr heat
+ !masout_ct - coefficient of turbulence (km) at cfl hight (m**2/s)
+ !masout_ckt - alft=kt/km ( kt-coefficient of turbulence for heat)
+
+ !> @brief Test - file selection for test
+
+ write(*,*) 'running code'
+
+ if (TEST==1) then
+ filename_in='MOSAiC.txt'
+ filename_out='out_MOSAiC.txt'
+ filename_in2='MOSAiC_zh.txt'
+ elseif (TEST==2) then
+ filename_in='Irgason1.txt'
+ filename_out='out_IRGASON1.txt'
+ filename_in2='IRGASON_zh.txt'
+ endif
+
+ open (1, file= filename_in, status ='old')
+ open (2, file=filename_out)
+ numst=0
+ do WHILE (ioer.eq.0)
+ read (1,*, iostat=ioer) data_in1%ws, data_in1%dt, data_in1%st, data_in1%dq
+ numst=numst+1
+ enddo
+ close (1)
+ numst=numst-1
+
+
+ allocate(data_inMAS%mas_w(numst))
+ allocate(data_inMAS%mas_dt(numst))
+ allocate(data_inMAS%mas_st(numst))
+ allocate(data_inMAS%mas_dq(numst))
+ allocate(data_inMAS%mas_cflh(numst))
+ allocate(data_inMAS%mas_z0in(numst))
+
+
+
+ allocate(data_outMAS%masout_zl(numst))
+ allocate(data_outMAS%masout_ri(numst))
+ allocate(data_outMAS%masout_re(numst))
+ allocate(data_outMAS%masout_lnzuzt(numst))
+ allocate(data_outMAS%masout_zu(numst))
+ allocate(data_outMAS%masout_ztout(numst))
+ allocate(data_outMAS%masout_rith(numst))
+ allocate(data_outMAS%masout_cm(numst))
+ allocate(data_outMAS%masout_ch(numst))
+ allocate(data_outMAS%masout_ct(numst))
+ allocate(data_outMAS%masout_ckt(numst))
+
+ open (11, file=filename_in2, status ='old')
+ open (1, file= filename_in, status ='old')
+ read (11, *) cflh, z0in
+ do i=1,numst
+ read (1,*) data_in1%ws, data_in1%dt, data_in1%st, data_in1%dq
+ data_inMAS%mas_w(i)=data_inMAS%mas_w(i)+data_in1%ws
+ data_inMAS%mas_dt(i)=data_inMAS%mas_dt(i)+data_in1%dt
+ data_inMAS%mas_st(i)=data_inMAS%mas_st(i)+data_in1%st
+ data_inMAS%mas_dq(i)=data_inMAS%mas_dq(i)+data_in1%dq
+ data_inMAS%mas_cflh(i)=cflh
+ data_inMAS%mas_z0in(i)=z0in
+ enddo
+
+
+ CALL surf_fluxMAS(data_inMAS%mas_w, data_inMAS%mas_dt, data_inMAS%mas_st, data_inMAS%mas_dq,&
+ data_inMAS%mas_cflh, data_inMAS%mas_z0in,&
+ data_outMAS%masout_zl, data_outMAS%masout_ri, data_outMAS%masout_re, data_outMAS%masout_lnzuzt,&
+ data_outMAS%masout_zu,data_outMAS%masout_ztout,data_outMAS%masout_rith,data_outMAS%masout_cm,&
+ data_outMAS%masout_ch,data_outMAS%masout_ct,data_outMAS%masout_ckt,&
+ data_par1, data_lutyp1,numst)
+
+ do i=1,numst
+ write (2,20) data_outMAS%masout_zl(i), data_outMAS%masout_ri(i), data_outMAS%masout_re(i), data_outMAS%masout_lnzuzt(i),&
+ data_outMAS%masout_zu(i), data_outMAS%masout_ztout(i), data_outMAS%masout_rith(i), data_outMAS%masout_cm(i),&
+ data_outMAS%masout_ch(i), data_outMAS%masout_ct(i), data_outMAS%masout_ckt(i)
+ enddo
+
+
+ deallocate(data_inMAS%mas_w)
+ deallocate(data_inMAS%mas_dt)
+ deallocate(data_inMAS%mas_st)
+ deallocate(data_inMAS%mas_dq)
+ deallocate(data_inMAS%mas_cflh)
+ deallocate(data_inMAS%mas_z0in)
- enddo
+ deallocate(data_outMAS%masout_zl)
+ deallocate(data_outMAS%masout_ri)
+ deallocate(data_outMAS%masout_re)
+ deallocate(data_outMAS%masout_lnzuzt)
+ deallocate(data_outMAS%masout_zu)
+ deallocate(data_outMAS%masout_ztout)
+ deallocate(data_outMAS%masout_rith)
+ deallocate(data_outMAS%masout_cm)
+ deallocate(data_outMAS%masout_ch)
+ deallocate(data_outMAS%masout_ct)
+ deallocate(data_outMAS%masout_ckt)
-100 continue
-10 format (4i4,5f7.1,f7.4,f7.1)
-20 format (4i4,5f7.1,f7.4,f7.1)
+ 10 format (f8.4,2x,f8.4)
+ 20 format (11(f10.4,3x))
stop
END PROGRAM
\ No newline at end of file
diff --git a/makefile b/makefile
index 1f03324822f3d197c8e05d77d3782636190e25d1..ef11984d0d2f2b5d67e8d3160ab6a97d97d02428 100644
--- a/makefile
+++ b/makefile
@@ -1,5 +1,5 @@
RUN = drag.exe
-COMPILER ?= intel
+COMPILER ?= gnu
FC_KEYS ?=
# set compiler
diff --git a/param.f90 b/param.f90
index ee8540e06daad09c4ec60b9796abfc5e99ef4176..467f4dec9963b55d45449ed1917fc911730a7407 100644
--- a/param.f90
+++ b/param.f90
@@ -1,27 +1,27 @@
-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
- 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 :: G0=1.2
- real, parameter :: AL2=(.14E0*(30.0E0**AN2))*(P4**.8E0)
- real, parameter :: A2=ALOG(H1*(G/ALFAM))
- real, parameter :: A3=BETAM*AN*A2
- real, parameter :: R0=.9E0/B4
+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
+ 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 :: g0=1.2
+ real, parameter :: al2=(.14e0*(30.0e0**an2))*(p4**.8e0)
+ real, parameter :: a2=alog(h1*(g/alfam))
+ real, parameter :: a3=betam*an*a2
+ real, parameter :: r0=.9e0/b4
!C* real, parameter :: AN5=(A6/A0)**4
!C* real, parameter :: D1=(2.0E0*G10-AN5*G4-SQRT((AN5*G4)**2+4.0E0*AN5*G10*(G10-G4)))/(2.0E0*G10**2)
!C* real, parameter :: Y10=(1.0E0-G4*D1)**.25E0
diff --git a/prmt.mod b/prmt.mod
deleted file mode 100644
index fbc5c106bbef1e8d13147a0421ec45fb0d9c8dd0..0000000000000000000000000000000000000000
Binary files a/prmt.mod and /dev/null differ