Commit 5ed15276 authored by ldemaine's avatar ldemaine
Browse files

Ajout de l'itération temporelle pour la convergence

parent 497f812a
......@@ -7983,484 +7983,6 @@ c$$$ close(1)
end
c$$$! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c$$$ subroutine slope_wind_backup(dset,xlon,xlat,zsurface,
c$$$ & zradius,zsurf_0,
c$$$ & temp,temp_gcm,tsurf,
c$$$ & Pl,sigma,ps_hr,
c$$$ & levhi,levlow,sheight,
c$$$ & utime,itimint,wl,wh,
c$$$ & zonwind,merwind,old_zonwind,old_merwind,
c$$$ & upslope_wind,crossslope_wind)
c$$$
c$$$
c$$$ implicit none
c$$$ include "constants_mcd.inc"
c$$$
c$$$! TODO: profil de la température
c$$$! TODO: calcul de K pour chaque couche
c$$$
c$$$
c$$$
c$$$
c$$$
c$$$! Inputs
c$$$ character*(*),intent(in) :: dset ! Path to MCD datafiles
c$$$
c$$$ real,intent(in) :: xlon
c$$$ real,intent(in) :: xlat
c$$$ real,intent(in) :: zsurface
c$$$ real,intent(in) :: zradius ! to compute get_slopes
c$$$ real,intent(in) :: zsurf_0
c$$$ real,intent(in) :: temp
c$$$ real,intent(in) :: Pl
c$$$ real,intent(in) :: temp_gcm(dimlevs)
c$$$ real,intent(in) :: tsurf
c$$$ real,intent(in) :: sigma(dimlevs),ps_hr
c$$$c$$$ !inutile pour le moment car temp déjà un imput
c$$$c$$$ real,intent(in) :: levweight(2)
c$$$ integer,intent(in) :: levhi
c$$$ integer,intent(in) :: levlow
c$$$c$$$ real,intent(in) :: pratio
c$$$ real,intent(in) :: sheight(dimlevs) ! altitude above surface of sigma levels
c$$$
c$$$ real,intent(in) :: utime !Universal time (0. to 24. hrs)=local time at lon=0
c$$$ integer,intent(in) :: itimint !seasonal interpolation flags
c$$$ real,intent(in) :: wl,wh !seasonal interpolation weights
c$$$
c$$$ real,intent(inout) :: zonwind
c$$$ real,intent(inout) :: merwind
c$$$
c$$$ real,intent(out) :: old_zonwind
c$$$ real,intent(out) :: old_merwind
c$$$
c$$$ real,intent(out) :: upslope_wind
c$$$ real,intent(out) :: crossslope_wind
c$$$ integer :: ier ! need intent(out)
c$$$
c$$$ real :: p_gcm(dimlevs)
c$$$
c$$$ real :: theta_s
c$$$ real :: psy_s
c$$$ real :: slopes_scale = 1
c$$$c$$$
c$$$c$$$ ! varibles pour la référence
c$$$ real :: zref
c$$$ integer :: levhi_0
c$$$ integer :: levlow_0
c$$$ real :: levweight_0
c$$$ real :: pratio_0
c$$$ integer :: l
c$$$ real temp_0
c$$$
c$$$
c$$$ real u_gcm(dimlevs),v_gcm(dimlevs)
c$$$ character*16 tmpname
c$$$ integer :: ierr
c$$$
c$$$ real :: Cd
c$$$ real :: Flotation
c$$$ real :: pKv(dimlevs),pKv_loc
c$$$ real Flotabis,Sciz
c$$$ real upwind_gcm(dimlevs),crosswind_gcm(dimlevs)
c$$$ integer it
c$$$c$$$ real :: UU
c$$$c$$$ integer :: i
c$$$c$$$ real :: uxu
c$$$c$$$ real :: uyu
c$$$c$$$
c$$$
c$$$
c$$$ real :: karman,g_0
c$$$ parameter (karman=0.4)
c$$$ parameter (g_0=3.72)
c$$$
c$$$ real pi,degtorad
c$$$ parameter (pi=3.14159265358979d0)
c$$$ parameter (degtorad=pi/180.0d0)
c$$$
c$$$ p_gcm = ps_hr*sigma
c$$$
c$$$ write(*,*) ' '
c$$$ write(*,*) '|--------------------------------------------------|'
c$$$ write(*,*) '| Input slope_wind |'
c$$$ write(*,*) '|--------------------------------------------------|'
c$$$ write(*,*) '| xlon, xlat:', xlon, xlat
c$$$ write(*,*) '| zsurface, zradius, zsurf_O:',zsurface,zradius,
c$$$ & zsurf_0
c$$$ write(*,*) '| levhi, levlow:', levhi, levlow
c$$$ write(*,*) '| utime:', utime
c$$$ write(*,*) '| itimint, wl, wh:', itimint, wl, wh
c$$$ write(*,*) '| zonwind, merwind:', zonwind, merwind
c$$$ write(*,*) '| temp, tsurf:', temp, tsurf
c$$$ write(*,*) '| temp_gcm:', temp_gcm
c$$$ write(*,*) '| sheight:', sheight
c$$$ write(*,*) '| sigma', sigma
c$$$ write(*,*) '| p_gcm', p_gcm
c$$$ write(*,*) '|--------------------------------------------------|'
c$$$ write(*,*) '| '
c$$$
c$$$c$$$ ! TODO: cas où l'altitude est inférieure à la première maille de mcd
c$$$
c$$$ upslope_wind=0
c$$$ crossslope_wind=0
c$$$
c$$$ old_merwind=merwind
c$$$ old_zonwind = zonwind
c$$$
c$$$ write(*,*) '| upslope_wind, crossslope_wind:'
c$$$ & ,upslope_wind, crossslope_wind
c$$$ write(*,*) '| old_merwind, merwind, old_zonwind, zonwind:'
c$$$ & ,old_merwind, merwind, old_zonwind, zonwind
c$$$ write(*,*) '| '
c$$$
c$$$ slopes_scale = 1
c$$$ call get_slopes(dset,xlon,xlat,slopes_scale,zradius,
c$$$ & theta_s,psy_s,ier)
c$$$
c$$$ zref = 10000*sin(theta_s*degtorad)
c$$$
c$$$ write(*,*) '| slopes_scale, theta_s, psy_s:',
c$$$ & slopes_scale, theta_s, psy_s
c$$$ write(*,*) '| zref:', zref
c$$$ write(*,*) '| '
c$$$
c$$$ if (zsurface.le.zref.and.zsurface.gt.0) then
c$$$ write(*,*) '| | if (zsurface.le.zref.and.zsurface.gt.0)'
c$$$ write(*,*) '|-| - - - - - - - - - - - - - - - - - - - - - - -'
c$$$ write(*,*) '| |'
c$$$
c$$$
c$$$
c$$$!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c$$$! calcul des valeurs pour la ref (sur le modèle de getsi)
c$$$!
c$$$! ATENTION pratio_0 non calculé dans les cas extrèmes pour l'instant
c$$$! (et non pris en compte dans les calculs de temp_0
c$$$
c$$$ if (zref.lt.sheight(1)) then
c$$$ write(*,*) '| | | if (zref.lt.sheight(1))'
c$$$ write(*,*) '|-|-|- - - - - - - - - - - - - - - - - - - - -'
c$$$! ! below the lowest layer
c$$$!c$$$ write(*,*) 'ATTENTION! zref trop bas'
c$$$ levhi_0=1
c$$$ levlow_0=1
c$$$ levweight_0=1.
c$$$ pratio_0 = 1.
c$$$
c$$$ elseif (zref.ge.sheight(dimlevs)) then
c$$$ write(*,*) '| | | elseif (zref.ge.sheight(dimlevs))'
c$$$ write(*,*) '|-|-|- - - - - - - - - - - - - - - - - - - - -'
c$$$!c$$$ write(*,*) 'ATTENTION! zref trop haut'
c$$$! ! above the top layer
c$$$ levhi_0=dimlevs
c$$$ levlow_0=dimlevs
c$$$ levweight_0=0.
c$$$ pratio_0 = 1.
c$$$ else
c$$$ write(*,*) '| | | else'
c$$$ write(*,*) '|-|-|- - - - - - - - - - - - - - - - - - - - -'
c$$$! ! general case: sheight(1)<=absheight<=sheight(dimlevs)
c$$$ do l=1,dimlevs-1
c$$$ if ((zref.ge.sheight(l))
c$$$ & .and.(zref.lt.sheight(l+1))) then
c$$$ levhi_0=l+1
c$$$ levlow_0=l
c$$$ levweight_0=(zref-sheight(levlow_0))
c$$$ & /(sheight(levhi_0)-sheight(levlow_0))
c$$$ pratio_0=1.
c$$$ endif
c$$$ enddo
c$$$
c$$$ endif
c$$$
c$$$ temp_0 = temp_gcm(levlow_0)
c$$$ & + (temp_gcm(levhi_0)-temp_gcm(levlow_0))*levweight_0
c$$$
c$$$ write(*,*) '| | | levhi_0, levlow_0, levweight_0, pratio_0:',
c$$$ & levhi_0,levlow_0,levweight_0,pratio_0
c$$$ write(*,*) '| | temp_0:', temp_0
c$$$ write(*,*) '| |'
c$$$
c$$$
c$$$
c$$$!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c$$$!
c$$$
c$$$ if (levhi.eq.1) then
c$$$ write(*,*) '| | | if (levhi.eq.1)'
c$$$ write(*,*) '| | |- - - - - - - - - - - - - - - - - - - - -'
c$$$ Flotation =g_0*(temp-temp_0)*sin(theta_s*degtorad)/temp_0
c$$$ Cd = (karman/log(zsurface/zsurf_0))**2
c$$$ upslope_wind = sign(sqrt(abs(Flotation / Cd)), Flotation)
c$$$ write(*,*) '| | | Flotation, Cd, upslope_wind:',
c$$$ & Flotation, Cd, upslope_wind
c$$$ write(*,*) '| | |'
c$$$
c$$$ else !(levhi.eq.1)
c$$$ write(*,*) '| | | if not (levhi.eq.1)'
c$$$ write(*,*) '| | |- - - - - - - - - - - - - - - - - - - - -'
c$$$! Calcul des valeurs au premier niveau gcm avec le coef Cd
c$$$ Flotation=g_0*(temp_gcm(1) - temp_0)*sin(theta_s*degtorad)
c$$$ & /temp_0
c$$$
c$$$ Cd = (karman/log(sheight(1)/zsurf_0))**2
c$$$
c$$$! Calcul du premier niveau de vent
c$$$ upwind_gcm(1)=sign(sqrt(abs(Flotation / Cd)), Flotation)
c$$$
c$$$ write(*,*) '| | | Flotation, Cd, upwind_gcm(1):',
c$$$ & Flotation, Cd, upwind_gcm(1)
c$$$ write(*,*) '| | |'
c$$$
c$$$! get wind profile at same location
c$$$ tmpname="u"
c$$$ call profi(u_gcm,xlon,xlat,utime,tmpname,ierr,itimint,
c$$$ & wl,wh,0,0,0)
c$$$ tmpname="v"
c$$$ call profi(v_gcm,xlon,xlat,utime,tmpname,ierr,itimint,
c$$$ & wl,wh,0,0,0)
c$$$ write(*,*) '| | | u_gcm:', u_gcm
c$$$ write(*,*) '| | |'
c$$$ write(*,*) '| | | v_gcm:', v_gcm
c$$$ write(*,*) '| | |'
c$$$
c$$$! Calcul de la première valeur de K (K(1))
c$$$ call vdif_k(sheight(1),0.,U_gcm(1),0.,V_gcm(1),0.,
c$$$ & Temp_gcm(1),tsurf,p_gcm(1),ps_hr,ps_hr,
c$$$ & zsurf_0,pKv(1))
c$$$ write(*,*) '| | | pKv(1):', pKv(1)
c$$$ write(*,*) '| | |'
c$$$
c$$$ if (levlow.eq.1) then
c$$$ write(*,*) '| | | | if levlow.eq.1'
c$$$ write(*,*) '|-|-|-| - - - - - - - - - - - - - - - - - -'
c$$$! K entre la dernière couche et la couche actuelle (pKv_loc)
c$$$ call vdif_k(sheight(2), sheight(1),
c$$$ & U_gcm(2), U_gcm(1),
c$$$ & V_gcm(2), V_gcm(1),
c$$$ & Temp_gcm(2), Temp_gcm(1),
c$$$ & p_gcm(2), p_gcm(1), ps_hr,
c$$$ & zsurf_0,
c$$$ & pKv_loc)
c$$$
c$$$ Flotation =g_0*(temp-temp_0)*sin(theta_s*degtorad)
c$$$ & /temp_0
c$$$
c$$$ Flotabis = (zsurface - sheight(1))*(zsurface)*Flotation
c$$$ & /(2*pKv_loc)
c$$$
c$$$ Sciz=(pkv(1)/pkv_loc)
c$$$ & *((zsurface-sheight(1))/(sheight(1)))
c$$$ & *(upwind_gcm(1))
c$$$
c$$$ upslope_wind = upwind_gcm(levlow) + Flotabis + Sciz
c$$$
c$$$ write(*,*) '| | | | pKv_loc', pKv_loc
c$$$ write(*,*) '| | | | Flotation, Flotabis, Sciz:',
c$$$ & Flotation, Flotabis, Sciz
c$$$ write(*,*) '| | | | upslope_wind:', upslope_wind
c$$$ write(*,*) '| | | | '
c$$$
c$$$ else !(levlow.eq.1)
c$$$ write(*,*) '| | | | if not levlow.eq.1'
c$$$ write(*,*) '|-|-|-| - - - - - - - - - - - - - - - - - -'
c$$$ do it=2,levlow
c$$$! K(it)
c$$$ write(*,*) '| | | | | do it =2,levlow'
c$$$ write(*,*) '|-|-|-|-|- - - - - - - - - - - - - - - -'
c$$$ call vdif_k(sheight(it), sheight(it-1),
c$$$ & U_gcm(it), U_gcm(it-1),
c$$$ & V_gcm(it),V_gcm(it-1),
c$$$ & Temp_gcm(it), Temp_gcm(it-1),
c$$$ & p_gcm(it), p_gcm(it-1), ps_hr,
c$$$ & zsurf_0,pKv(it))
c$$$
c$$$ Flotation = g_0*(temp_gcm(it)-temp_0)
c$$$ & *sin(theta_s*degtorad)/temp_0
c$$$
c$$$ if (it.eq.2) then
c$$$ Sciz = (pkv(1)/pkv(2))
c$$$ & *((sheight(2)-sheight(1))/sheight(1))
c$$$ & *upwind_gcm(1)
c$$$
c$$$ Flotabis = (sheight(2)-sheight(1))*sheight(2)
c$$$ & *Flotation/(2*pkv(2))
c$$$
c$$$ else !(it.eq.2)
c$$$ Sciz = (pkv(it-1)/pkv(it))
c$$$ & *((sheight(it)-sheight(it-1))
c$$$ & /(sheight(it-1)-sheight(it-2)))
c$$$ & *(upwind_gcm(it-1)-upwind_gcm(it-2))
c$$$
c$$$ Flotabis = (sheight(it)-sheight(it-1))
c$$$ & *(sheight(it)-sheight(it-2))
c$$$ & *Flotation/(2*pkv(it))
c$$$
c$$$ endif !(it.eq.2)
c$$$
c$$$ upwind_gcm(it) = upwind_gcm(it-1) + Flotabis + Sciz
c$$$
c$$$ write(*,*) '| | | | | pKv(it):',pKv(it)
c$$$ write(*,*) '| | | | | Flotation, Flotabis, Sciz:',
c$$$ & Flotation, Flotabis, Sciz
c$$$ write(*,*) '| | | | | upwind_gcm(it):',upwind_gcm(it)
c$$$ write(*,*) '| | | | |'
c$$$
c$$$ enddo !it=2,levlow
c$$$
c$$$! K entre la dernière couche et la couche actuelle (pKv_loc)
c$$$ call vdif_k(zsurface, sheight(levlow),
c$$$ & merwind, U_gcm(levlow),
c$$$ & zonwind, V_gcm(levlow),
c$$$ & temp, Temp_gcm(levlow),
c$$$ & Pl, p_gcm(levlow), ps_hr,
c$$$ & zsurf_0,
c$$$ & pKv_loc)
c$$$
c$$$ Flotation =g_0*(temp-temp_0)*sin(theta_s*degtorad)
c$$$ & /temp_0
c$$$
c$$$ Flotabis = (zsurface - sheight(levlow))
c$$$ & *(zsurface - sheight(levlow-1))
c$$$ & *Flotation/(2*pKv_loc)
c$$$
c$$$ Sciz = (pkv(levlow)/pkv_loc)
c$$$ & *((zsurface-sheight(levlow))
c$$$ & /(sheight(levlow)-sheight(levlow-1)))
c$$$ & *(upwind_gcm(levlow)-upwind_gcm(levlow-1))
c$$$
c$$$ upslope_wind = upwind_gcm(levlow) + Flotabis + Sciz
c$$$
c$$$ write(*,*) '| | | | pKv_loc', pKv_loc
c$$$ write(*,*) '| | | | Flotation, Flotabis, Sciz:',
c$$$ & Flotation, Flotabis, Sciz
c$$$ write(*,*) '| | | | upslope_wind:', upslope_wind
c$$$ write(*,*) '| | | | '
c$$$
c$$$ endif !(levlow.eq.1)
c$$$ endif !(levhi.eq.1)
c$$$
c$$$ endif !(zsurface.le.zref.and.zsurface.gt.0)
c$$$
c$$$ merwind = old_merwind
c$$$ & + upslope_wind*cos(psy_s*degtorad)
c$$$ & + crossslope_wind*sin(psy_s*degtorad)
c$$$
c$$$ zonwind = old_zonwind
c$$$ & + upslope_wind*sin(psy_s*degtorad)
c$$$ & - crossslope_wind*cos(psy_s*degtorad)
c$$$
c$$$ write(*,*) '| merwind, zonwind:', merwind, zonwind
c$$$ return
c$$$
c$$$ end
c$$$
c$$$
c$$$! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c$$$ subroutine vdif_k_backup(z_hi,z_low,
c$$$ & U_hi,u_low,
c$$$ & V_hi,V_low,
c$$$ & temp_hi,temp_low,
c$$$ & P_hi,P_low,Psurf,
c$$$ & zsurf_0,
c$$$ & pKv)
c$$$
c$$$ implicit none
c$$$ include "constants_mcd.inc"
c$$$
c$$$ real,intent(in) :: z_hi,z_low
c$$$ real,intent(in) :: U_hi,U_low
c$$$ real,intent(in) :: V_hi,V_low
c$$$ real,intent(in) :: Temp_low,Temp_hi
c$$$ real,intent(in) :: P_hi,P_low
c$$$ real,intent(in) :: Psurf
c$$$ real,intent(in) :: zsurf_0
c$$$
c$$$ real,intent(out) :: pKv
c$$$
c$$$
c$$$ real :: zdu,zdv,zdz
c$$$ real :: zdvodz2
c$$$
c$$$ real :: rcp
c$$$ real :: T_pot_low,T_pot_hi
c$$$ real :: T_pot_mid,dT_pot_dz
c$$$
c$$$ real :: z1
c$$$ real :: Ri
c$$$
c$$$ real :: lmix,lmixmin
c$$$ real :: emin_turb
c$$$
c$$$ real :: karman, g_0
c$$$ parameter (karman=0.4)
c$$$ parameter (g_0=3.72)
c$$$
c$$$ write(*,*) ' > ********************************************'
c$$$ write(*,*) ' > Input vdif_k'
c$$$ write(*,*) ' > ********************************************'
c$$$ write(*,*) ' > z_0, z_hi, z_low:',
c$$$ & zsurf_0, z_hi, z_low
c$$$ write(*,*) ' > U_low, U_hi, V_low, V_hi:',
c$$$ & U_low, U_hi, V_low, V_hi
c$$$ write(*,*) ' > temp_low, Temp_hi:',
c$$$ & temp_low, Temp_hi
c$$$ write(*,*) ' > P_low, P_hi, Psurf:',
c$$$ & P_low, P_hi, Psurf
c$$$ write(*,*) ' >'
c$$$
c$$$ lmixmin = 30.
c$$$ emin_turb = 1.e-8
c$$$
c$$$
c$$$ z1 = (z_hi + z_low)/2 + zsurf_0
c$$$ lmix = karman*z1/(1.+karman*z1/lmixmin)
c$$$
c$$$ zdu = U_hi - U_low
c$$$ zdv = V_hi - V_low
c$$$ zdz = z_hi - z_low
c$$$ zdvodz2 = (zdu**2 + zdv**2)/(zdz**2)
c$$$ write(*,*) ' > z1, lmix, zdu, zdv, zdz, zdvodz2:',
c$$$ & z1, lmix, zdu, zdv, zdz, zdvodz2
c$$$ write(*,*) ' >'
c$$$
c$$$! need declaration
c$$$ rcp=0.259 ! on considère que gamma = 1.35 avec rcp = (1-gamma)/gamma, ce qui est valable à 3% près en dessous de 5000m
c$$$ T_pot_hi=Temp_hi*(Psurf/P_hi)**rcp
c$$$ T_pot_low=Temp_low*(Psurf/P_low)**rcp
c$$$ T_pot_mid = (T_pot_hi + T_pot_low)/2.
c$$$ dT_pot_dz = (T_pot_hi - T_pot_low)/zdz
c$$$ write(*,*) ' > T_pot_hi, T_pot_low, T_pot_mid, dT_pot_dz:',
c$$$ & T_pot_hi, T_pot_low, T_pot_mid, dT_pot_dz
c$$$ write(*,*) ' >'
c$$$
c$$$ write(*,*) ' >'
c$$$ write(*,*) ' > zdvodz2:', zdvodz2
c$$$ write(*,*) ' > lmix', lmix
c$$$ write(*,*) ' >'
c$$$ write(*,*) ' > dT_pot_dz:', dT_pot_dz
c$$$ write(*,*) ' > g_0:', g_0
c$$$ write(*,*) ' >'
c$$$
c$$$ if(zdvodz2.lt.1.e-5) then
c$$$ pKv = lmix*sqrt(emin_turb)
c$$$ else
c$$$ Ri = g_0 * dT_pot_dz /(T_pot_mid*zdvodz2)
c$$$
c$$$ write(*,*) ' > * Ri:', Ri
c$$$ write(*,*) ' > 1 - 2.5 Ri:', 1.-2.5*Ri
c$$$ write(*,*) ' >'
c$$$
c$$$ pKv = lmix*sqrt(max(lmix**2*zdvodz2*(1-2.5*Ri),emin_turb))
c$$$ endif
c$$$ write(*,*) ' > pKv:', pKv
c$$$ write(*,*) ' > ********************************************'
c$$$
c$$$ return
c$$$ end
c$$$
! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine slope_wind(dset,xlon,xlat,zsurface,
& zradius,zsurf_0,
......@@ -8475,11 +7997,6 @@ c$$$
implicit none
include "constants_mcd.inc"
! TODO: profil de la température
! TODO: calcul de K pour chaque couche
! Inputs
......@@ -8540,9 +8057,14 @@ c$$$ ! varibles pour la référence
real :: Cd
real :: Flotation
real :: pKv(dimlevs),pKv_loc
real Flotabis,Sciz
real upwind_gcm(dimlevs),crosswind_gcm(dimlevs)
integer it
real :: Flotabis
real :: Sciz_u, Sciz_v
real :: Corio_u, Corio_v
real :: upwind_gcm(dimlevs),crosswind_gcm(dimlevs)
real :: u_prev(dimlevs),v_prev(dimlevs)
integer :: it
real :: f
real :: dt
c$$$ real :: UU
c$$$ integer :: i
c$$$ real :: uxu
......@@ -8712,7 +8234,11 @@ c$$$ ! TODO: cas où l'altitude est inférieure à la première maille de mcd
!
do l=1,2
dt = 0.01
f = sin(degtorad*xlat)*4*pi/(24.6*3600)
do l=1,10
write(*,*) '| | | | (do l=1,2), l=',l
write(*,*) '|-|-|-|- - - - - - - - - - - - - - - - - -'
......@@ -8727,31 +8253,67 @@ c$$$ ! TODO: cas où l'altitude est inférieure à la première maille de mcd
write(*,*) '| | | | pKv', pKv
write(*,*) '| | | |'
c$$$ do it=2,levhi
c$$$ Flotation = g_0*(temp_gcm(it)-temp_0)
c$$$ & *sin(theta_s*degtorad)/temp_0
c$$$
c$$$ if (it.eq.2) then
c$$$ Sciz = (pkv(1)/pkv(2))
c$$$ & *((sheight(2)-sheight(1))/sheight(1))
c$$$ & *upwind_gcm(1)
c$$$
c$$$ Flotabis = (sheight(2)-sheight(1))*sheight(2)
c$$$ & *Flotation/(2*pkv(2))
c$$$
c$$$ else !(it.eq.2)
c$$$ Sciz = (pkv(it-1)/pkv(it))
c$$$ & *((sheight(it)-sheight(it-1))
c$$$ & /(sheight(it-1)-sheight(it-2)))
c$$$ & *(upwind_gcm(it-1)-upwind_gcm(it-2))
c$$$
c$$$ Flotabis = (sheight(it)-sheight(it-1))
c$$$ & *(sheight(it)-sheight(it-2))
c$$$ & *Flotation/(2*pkv(it))
c$$$
c$$$ endif !(it.eq.2)
c$$$
c$$$ upwind_gcm(it) = upwind_gcm(it-1) - Flotabis + Sciz
c$$$
c$$$ u_gcm(it) = old_u_gcm(it)
c$$$ & + upwind_gcm(it)*cos(psy_s*degtorad)
c$$$ & + crosswind_gcm(it)*sin(psy_s*degtorad)
c$$$
c$$$ v_gcm(it) = old_v_gcm(it)
c$$$ & + upwind_gcm(it)*cos(psy_s*degtorad)
c$$$ & - crosswind_gcm(it)*sin(psy_s*degtorad)
c$$$
c$$$ enddo
do it=2,levhi
Flotation = g_0*(temp_gcm(it)-temp_0)
& *sin(theta_s*degtorad)/temp_0
if (it.eq.2) then
Sciz = (pkv(1)/pkv(2))
& *((sheight(2)-sheight(1))/sheight(1))
& *upwind_gcm(1)
Sciz_u = (2/(sheight(it+1)-sheight(it-1)))
& *(pKv(it+1)*(u_prev(it+1)-u_prev(it))
& /(sheight(it+1)-sheight(it))
& - pKv(it)*(u_prev(it)-u_prev(it-1))
& /(sheight(it)-sheight(it-1)))