Skip to content
Snippets Groups Projects
Commit 7bd79379 authored by Evgeny Mortikov's avatar Evgeny Mortikov
Browse files

update to include MOSAiC data; array calculations

parent 593c78e4
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
.gitignore 0 → 100644
+ 1
− 0
View file @ 7bd79379
/drag.exe
compile2.sh 100644 → 100755
+ 0
− 1
View file @ 7bd79379
......@@ -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
......
drag3.f90
+ 105
− 40
View file @ 7bd79379
......@@ -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
......
inputdata.f90
+ 8
− 1
View file @ 7bd79379
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: =
......
main_drag.f90
+ 151
− 14
View file @ 7bd79379
......@@ -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
makefile
+ 1
− 1
View file @ 7bd79379
RUN = drag.exe
COMPILER ?= intel
COMPILER ?= gnu
FC_KEYS ?=
# set compiler
......
param.f90
+ 24
− 24
View file @ 7bd79379
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
......
prmt.mod deleted 100644 → 0
+ 0
− 0
View file @ 593c78e4
File deleted
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment