From 7bd7937964d6c0ffd9ef4e0c7dc99a67130b14f7 Mon Sep 17 00:00:00 2001
From: Evgeny Mortikov <evgeny.mortikov@gmail.com>
Date: Fri, 15 Dec 2023 21:44:22 +0300
Subject: [PATCH] update to include MOSAiC data; array calculations
---
.gitignore | 1 +
compile2.sh | 1 -
drag3.f90 | 145 ++++++++++++++++++++++++++++++++------------
inputdata.f90 | 9 ++-
main_drag.f90 | 165 +++++++++++++++++++++++++++++++++++++++++++++-----
makefile | 2 +-
param.f90 | 48 +++++++--------
prmt.mod | Bin 601 -> 0 bytes
8 files changed, 290 insertions(+), 81 deletions(-)
create mode 100644 .gitignore
mode change 100644 => 100755 compile2.sh
delete mode 100644 prmt.mod
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..b2e5ffe
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1 @@
+/drag.exe
diff --git a/compile2.sh b/compile2.sh
old mode 100644
new mode 100755
index e1433f8..fac0309
--- 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 039a999..41b5bd3 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 8fb6737..3543420 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 6c809af..6e4faff 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 1f03324..ef11984 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 ee8540e..467f4de 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
GIT binary patch
literal 0
HcmV?d00001
literal 601
zcmV-f0;c^RiwFP!000001J#&aZ<{a_hVT6=&Yj_+lWjf%mx#~`tuYZLn|f7cEs>%o
zq$0H1{{0TYhNLTO+IC|>#5ut>=i&I{H8<B;UgVL4{o(0#SHfFaAGe1JDEBBl*5$e>
zpI}oT_V7~gn@^jN5zDgSf1E@cHbL`)p8&dDQ=sjW@j`)vWi6vTS|ml1!%BXZ*;jd$
z$|8|P?=JH!)(%M*%iA<gi>tfnO9Bj7_sF)&4L0Gnf8n8hjOFCFW#UZU6%jUUk?G=;
z)01V{s-`lTO&p-lI|B&BID*?E?9FEIkm1-Ez`$pt3_j4^`nkkmTK-*+BFce(ARz<F
zZ%EiE$P7Wgy>3j%M#%X%M}Qn<)qFPq#Pxj<xZ^0HLNS3QbcE*z<18G@;2=Q+Yg;w-
zwmNPf?EwhmP3EgS`J<7duWcCDD;XDQCLzl6=>8#F;aho_C(A_66B*xkF`QSxqi^UO
z)dv1c%p7R7*_bdx2}g{-dR(xp-P!~UKCxPGr^DEDBS0pGg4=?WZw6pYrNiqej7Y~+
zHvE?3P#HkMr{>pLhw)-UI(RVAuL1=ki~yOKc5(TxYGFb?#5je(+qhmE05PEn1wTRj
zFo41thY)@A7Z8hVzPg<hPCdby4*F0S^~&l}7`>0A{RnXS1i(OaJ81@v+K)1D)vk_#
zcd<8O;Oi|0gPB(8gh4RKTn3?%iL8}lN`0ZtwS+-i3_1WkM=hq<kX~5T9Zyf}X#2zL
nh~8|NV$cipI7K@oU2jJ+ruz)qT)%U)g-w0}!l@Ipo)G{5Fu)>l
--
GitLab