********************************************************************* c Pogoreltsev November 2015 C Rf_c=0.4 C T_0, U_0, lower boundary fields from the UKMO data C C --- INPUT: nb, igit, zoben, rm(kgit), rhos (kg/m^3), skh C-------------------------------------------------------------------- C --- OUTPUT: C dz,defi,zoben,zunten,nstep,dt,dt2,cor(nb+1),tgfia(nb+1),... C dx(nb+1),phi(nb+1) C rou(kgit),row(kgit),alfa(kgit),drag(nb,kgit,1),drag(nb,kgit,2) C dichte(kgit), C a1h2o(igit,nb), C skhnuv(9),skhnir(4),e1nuv(9),e1nir(4) --> strobel.f ircool.f C------------------------------------------------------------------- subroutine init_60_field_sig c c-------------------------------------- INITIALIZATION FOR MUAM_60 c include 'com_main.fc' include 'data_TMP_SPW.f' include 'data_GH_SPW.f' include 'data_H2O_3d.f' c c -------------------------------------------------------------- c ntime anzahl der zeitschritte pro stunde (dt=3600/ntime c ntime ist so zu waehlen, dass dt ganz zahl) c ncom aktueller zeitschritt c nsec aktueller zeitpunkt c nend letzter prognostizierter zeitschritt c nout anzahl der zeitschritte bei denen jeweils auslagerung c von zwischenergebnissen auf band oder platte c nprint anzahl der zeitschritte bei denen jeweils ausdruck c von zwischenergebnissen erfolgt c igit anzahl der gitterpunkte laengs eines breitenkreises c kgit anzahl der modellflaechen c nwelle(1-3) gibt an, welche wellen im forcing enthalten sind c nt0 t0 in forcing function c nfase(1-3) phase der wellen 1-3 c jforc1 breitenkreis bei dem forcing beginnt c jforc2 breitenkreis bei dem forcing endet c nfinul(1-3) max. amplitude des forcing c rnn : vert. profile von n2, o!!!!!!!!!!!!!!!!!!!!!!!!!, o2, no c rm : vert. profil des molekulargewichts c rl (kgit+1) molekulare waermeleitfaehigkeit c rx (kgit+1) dynamische viskositaet c xmt (kgit+1) xmt=m'/tau' c rou (kgit) z-koordinate im niveau k c row (kgit) z-koordinate im niveau k-0.5 c---------------------------------------------------------------------- c REAL T0_NCEP(nb),G0_NCEP(nb) REAL*4 wk(nb,19) REAL gh1a(36),gh1f(36),gh2a(36),gh2f(36),gh3a(36),gh3f(36), $ gh4a(36),gh4f(36),gh5a(36),gh5f(36) REAL tp1a(36),tp1f(36),tp2a(36),tp2f(36),tp3a(36),tp3f(36), $ tp4a(36),tp4f(36),tp5a(36),tp5f(36) REAL dkli(9),xncir(4),rh2o(64,36) REAL alfaNC(60),sigma1(60),sigma2(60) c REAL*4 aus_0211(igit,nb,19,3) REAL*4 gh_ave(igit,nb),tp_ave(igit,nb) REAL*4 s1u_m1(5,nb,23,4),s1u_m2(5,nb,23,4), & s1u_m3(5,nb,23,4),s1u_m4(5,nb,23,4) REAL*4 s2_m1 (5,nb,23,4),s2_m2 (5,nb,23,4), & s2_m3 (5,nb,23,4),s2_m4 (5,nb,23,4) REAL*4 ht_m1 (3,nb,10,4),ht_m2 (3,nb,10,4), & ht_m3 (3,nb,10,4),ht_m4 (3,nb,10,4), & ht_m0 (2,nb,10,2) c------------------------------------------------ CONDUCTIVITIES REAL*4 s2_a(nb,23),s1u_a(nb,23),s1u_a0(nb,23),s2_a0(nb,23) REAL*4 ht_a(nb,10),ht_a0(nb,10) data a /6366197./ C-------------------------------------------------- K_i/D_i, Liou (1992) C---------------- pressure factors data dkli /0.226,0.25,0.54,0.52,0.43,0.62,0.62,0.51,0.87/ c------------------------------ coeff. For h2o column in ircool_newCO2.f c data xncir /0.5,0.5,0.5,0.5/ ! xnc data xncir /1.,1.,1.,1./ ! xnc c--------------------------------------------------- PHYSICAL PARAMETERS: c c---------- daily averaged - /1.5 c data sigma1/ 35*0.0, !64 levels c $0.02,0.024,0.08,0.29,0.48,1.22,4.58,11.2,23.0,38.2,47.9,54.8, c $13*60.5/ c---------- daily averaged - /1.5 c data sigma2/ 35*0.0, !64 level c $0.2,0.8,1.8,3.2,4.9,7.7,13.3,20.0,23.7,21.1,14.5,9.4, c $13*6.0/ c-------------------- Newtonian cooling prfile c Zhu, (1993) JAS, v. 50, No 17, p. 3008-3021. data alfaNC/3*0.070,0.08,0.110,0.124,0.176,0.277,0.390,0.533, $ 0.674,0.859,1.140,1.430,1.690,1.950,2.200,2.360, $ 2.360,2.350,2.390,2.520,2.650,2.750,2.760,2.760, $ 2.750,2.660,2.390,2.100,2.230,2.810,3.660,4.650, $ 5.180,6.700,7.310,7.210,22*6.890/ pi=2.*asin(1.) Open (55, file='ZW_Ely_Mar.dx', $ access='direct',status='unknown',recl=nb*17*1) Open (45, file='TP_Ely_Mar.dx', $ access='direct',status='unknown',recl=nb*17*1) c OPEN (202, FILE='UVT_0211.dx', c $ form='unformatted', c $ access='direct', status='unknown',recl=igit*nb*19*3) OPEN (203, FILE='GH_1000_Ely_Mar.dx', $ form='unformatted', $ access='direct', status='unknown',recl=igit*nb) OPEN (204, FILE='TP_1000_Ely_Mar.dx', $ form='unformatted', $ access='direct', status='unknown',recl=igit*nb) read(55,rec=1) U_UKMO read(45,rec=1) T_NCEP 211 format(1x,3e13.4) 1311 format(1x,F9.3,f10.3,F10.1) 212 format(1x,3e13.4) 213 format(1x,2e13.4) c-----------INPUT for the ionospheric conductivities ------------ c----------------- S2 -------------- open (11, FILE='s2_24_m1.dx' , $ form='unformatted', * access='direct', status='old',recl=5*nb*23*4) open (12, FILE='s2_24_m2.dx' , $ form='unformatted', * access='direct', status='old',recl=5*nb*23*4) c open (13, FILE='s2_24_m3.dx' , c $ form='unformatted', c * access='direct', status='old',recl=5*nb*23*4) c open (14, FILE='s2_24_m4.dx' , c $ form='unformatted', c * access='direct', status='old',recl=5*nb*23*4) c- -------------- S1U -------------- open (31, FILE='s1u_24_m1.dx' , $ form='unformatted', * access='direct', status='old',recl=5*nb*23*4) open (32, FILE='s1u_24_m2.dx' , $ form='unformatted', * access='direct', status='old',recl=5*nb*23*4) c open (33, FILE='s1u_24_m3.dx' , c $ form='unformatted', c * access='direct', status='old',recl=5*nb*23*4) c open (34, FILE='s1u_24_m4.dx' , c $ form='unformatted', c * access='direct', status='old',recl=5*nb*23*4) c-------------------------------------------------------- open (40, FILE='LH_m0_Mar.dx' , $ form='unformatted', * access='direct', status='old',recl=2*nb*10*2) open (41, FILE='lh_m1_av_Mar.dx' , $ form='unformatted', * access='direct', status='old',recl=3*nb*10*4) open (42, FILE='lh_m2_av_Mar.dx' , $ form='unformatted', * access='direct', status='old',recl=3*nb*10*4) open (43, FILE='lh_m3_av_Mar.dx' , $ form='unformatted', * access='direct', status='old',recl=3*nb*10*4) open (44, FILE='lh_m4_av_Mar.dx' , $ form='unformatted', * access='direct', status='old',recl=3*nb*10*4) c-------------------- read input data on harmonics READ(11,REC=1) s2_m1 READ(12,REC=1) s2_m2 c READ(13,REC=1) s2_m3 c READ(14,REC=1) s2_m4 READ(31,REC=1) s1u_m1 READ(32,REC=1) s1u_m2 c READ(33,REC=1) s1u_m3 c READ(34,REC=1) s1u_m4 READ(40,REC=1) ht_m0 READ(41,REC=1) ht_m1 READ(42,REC=1) ht_m2 READ(43,REC=1) ht_m3 READ(44,REC=1) ht_m4 c-------------------------------- Amplitudes and Phases do k=1,23 do j=1,nb c------------------------------------- As1m1w24(j,k)=s1u_m1(1,j,k,1) Fs1m1w24(j,k)=s1u_m1(1,j,k,3) As1m1w12(j,k)=s1u_m1(2,j,k,1) Fs1m1w12(j,k)=s1u_m1(2,j,k,3) As1m1w08(j,k)=s1u_m1(3,j,k,1) Fs1m1w08(j,k)=s1u_m1(3,j,k,3) As1m1w06(j,k)=s1u_m1(4,j,k,1) Fs1m1w06(j,k)=s1u_m1(4,j,k,3) As1m1spw(j,k)=s1u_m1(5,j,k,1) Fs1m1spw(j,k)=s1u_m1(5,j,k,3) c-------------------------------------- As1m1e24(j,k)=s1u_m1(1,j,k,2) Fs1m1e24(j,k)=s1u_m1(1,j,k,4) As1m1e12(j,k)=s1u_m1(2,j,k,2) Fs1m1e12(j,k)=s1u_m1(2,j,k,4) As1m1e08(j,k)=s1u_m1(3,j,k,2) Fs1m1e08(j,k)=s1u_m1(3,j,k,4) As1m1e06(j,k)=s1u_m1(4,j,k,2) Fs1m1e06(j,k)=s1u_m1(4,j,k,4) c-------------------------------------- As1m2w24(j,k)=s1u_m2(1,j,k,1) Fs1m2w24(j,k)=s1u_m2(1,j,k,3) As1m2w12(j,k)=s1u_m2(2,j,k,1) Fs1m2w12(j,k)=s1u_m2(2,j,k,3) As1m2w08(j,k)=s1u_m2(3,j,k,1) Fs1m2w08(j,k)=s1u_m2(3,j,k,3) As1m2w06(j,k)=s1u_m2(4,j,k,1) Fs1m2w06(j,k)=s1u_m2(4,j,k,3) As1m2spw(j,k)=s1u_m2(5,j,k,1) Fs1m2spw(j,k)=s1u_m2(5,j,k,3) c-------------------------------------- As1m1e24(j,k)=s1u_m2(1,j,k,2) Fs1m1e24(j,k)=s1u_m2(1,j,k,4) As1m1e12(j,k)=s1u_m2(2,j,k,2) Fs1m1e12(j,k)=s1u_m2(2,j,k,4) As1m1e08(j,k)=s1u_m2(3,j,k,2) Fs1m1e08(j,k)=s1u_m2(3,j,k,4) As1m1e06(j,k)=s1u_m2(4,j,k,2) Fs1m1e06(j,k)=s1u_m2(4,j,k,4) c###################################### As2m1w24(j,k)=s2_m1(1,j,k,1) Fs2m1w24(j,k)=s2_m1(1,j,k,3) As2m1w12(j,k)=s2_m1(2,j,k,1) Fs2m1w12(j,k)=s2_m1(2,j,k,3) As2m1w08(j,k)=s2_m1(3,j,k,1) Fs2m1w08(j,k)=s2_m1(3,j,k,3) As2m1w06(j,k)=s2_m1(4,j,k,1) Fs2m1w06(j,k)=s2_m1(4,j,k,3) As2m1spw(j,k)=s2_m1(5,j,k,1) Fs2m1spw(j,k)=s2_m1(5,j,k,3) c-------------------------------------- As2m1e24(j,k)=s2_m1(1,j,k,2) Fs2m1e24(j,k)=s2_m1(1,j,k,4) As2m1e12(j,k)=s2_m1(2,j,k,2) Fs2m1e12(j,k)=s2_m1(2,j,k,4) As2m1e08(j,k)=s2_m1(3,j,k,2) Fs2m1e08(j,k)=s2_m1(3,j,k,4) As2m1e06(j,k)=s2_m1(4,j,k,2) Fs2m1e06(j,k)=s2_m1(4,j,k,4) c-------------------------------------- As2m2w24(j,k)=s2_m2(1,j,k,1) Fs2m2w24(j,k)=s2_m2(1,j,k,3) As2m2w12(j,k)=s2_m2(2,j,k,1) Fs2m2w12(j,k)=s2_m2(2,j,k,3) As2m2w08(j,k)=s2_m2(3,j,k,1) Fs2m2w08(j,k)=s2_m2(3,j,k,3) As2m2w06(j,k)=s2_m2(4,j,k,1) Fs2m2w06(j,k)=s2_m2(4,j,k,3) As2m2spw(j,k)=s2_m2(5,j,k,1) Fs2m2spw(j,k)=s2_m2(5,j,k,3) c-------------------------------------- As2m2e24(j,k)=s2_m2(1,j,k,2) Fs2m2e24(j,k)=s2_m2(1,j,k,4) As2m2e12(j,k)=s2_m2(2,j,k,2) Fs2m2e12(j,k)=s2_m2(2,j,k,4) As2m2e08(j,k)=s2_m2(3,j,k,2) Fs2m2e08(j,k)=s2_m2(3,j,k,4) As2m2e06(j,k)=s2_m2(4,j,k,2) Fs2m2e06(j,k)=s2_m2(4,j,k,4) c-------------------------------------- end do end do c--------------------------------------------------- do k=1,10 do j=1,nb c------------------------------------- Alhm0_24(j,k)=ht_m0(1,j,k,1) Flhm0_24(j,k)=ht_m0(1,j,k,2) Alhm0_12(j,k)=ht_m0(2,j,k,1) Flhm0_12(j,k)=ht_m0(2,j,k,2) c------------------------------------- Alhm1w24(j,k)=ht_m1(1,j,k,1) Flhm1w24(j,k)=ht_m1(1,j,k,3) Alhm1w12(j,k)=ht_m1(2,j,k,1) Flhm1w12(j,k)=ht_m1(2,j,k,3) Alhm1e24(j,k)=ht_m1(1,j,k,2) Flhm1e24(j,k)=ht_m1(1,j,k,4) Alhm1e12(j,k)=ht_m1(2,j,k,2) Flhm1e12(j,k)=ht_m1(2,j,k,4) c Alhm1spw(j,k)=ht_m1(5,j,k,1) c Flhm1spw(j,k)=ht_m1(5,j,k,3) c-------------------------------------- Alhm2w24(j,k)=ht_m2(1,j,k,1) Flhm2w24(j,k)=ht_m2(1,j,k,3) Alhm2w12(j,k)=ht_m2(2,j,k,1) Flhm2w12(j,k)=ht_m2(2,j,k,3) c Alhm2spw(j,k)=ht_m2(5,j,k,1) c Flhm2spw(j,k)=ht_m2(5,j,k,3) c-------------------------------------- Alhm1e24(j,k)=ht_m2(1,j,k,2) Flhm1e24(j,k)=ht_m2(1,j,k,4) Alhm1e12(j,k)=ht_m2(2,j,k,2) Flhm1e12(j,k)=ht_m2(2,j,k,4) c------------------------------------- Alhm3w24(j,k)=ht_m3(1,j,k,1) Flhm3w24(j,k)=ht_m3(1,j,k,3) Alhm3w12(j,k)=ht_m3(2,j,k,1) Flhm3w12(j,k)=ht_m3(2,j,k,3) Alhm3e24(j,k)=ht_m3(1,j,k,2) Flhm3e24(j,k)=ht_m3(1,j,k,4) Alhm3e12(j,k)=ht_m3(2,j,k,2) Flhm3e12(j,k)=ht_m3(2,j,k,4) c------------------------------------- Alhm4w24(j,k)=ht_m4(1,j,k,1) Flhm4w24(j,k)=ht_m4(1,j,k,3) Alhm4w12(j,k)=ht_m4(2,j,k,1) Flhm4w12(j,k)=ht_m4(2,j,k,3) Alhm4e24(j,k)=ht_m4(1,j,k,2) Flhm4e24(j,k)=ht_m4(1,j,k,4) Alhm4e12(j,k)=ht_m4(2,j,k,2) Flhm4e12(j,k)=ht_m4(2,j,k,4) c###################################### end do end do c---------------------------- read T for March-April c do l=91,150 c read(202,rec=l) aus_0211 c do k=1,19 c do j=1,nb c do i=1,igit c T_UKMO_TIME(i,j,k,l-90)=aus_0211(i,j,k,3) c end do c end do c end do c end do c---------------------------------------- READ LOWER BOUNDARY CONDITIONS read(203,rec=1)gh_ave read(204,rec=1)tp_ave 2010 format(64F7.1) day=65. td1=2.*pi*day/365. td2=2.*td1 td3=3.*td1 td4=4.*td1 do 200 j=1,nb c-------------------------------------------------------------------- c rh2o(j)=h2oav_00(j) + h2oav_s1(j)*SIN(td1)+ h2oav_c1(j)*COS(td1)+ c $ h2oav_s2(j)*SIN(td2)+ h2oav_c2(j)*COS(td2)+ c $ h2oav_s3(j)*SIN(td3)+ h2oav_c3(j)*COS(td3)+ c $ h2oav_s4(j)*SIN(td4)+ h2oav_c4(j)*COS(td4) c-------------------------------------------------------------------- c H_wv(j)=hh2oav_00(j)+hh2oav_s1(j)*SIN(td1)+hh2oav_c1(j)*COS(td1)+ c $ hh2oav_s2(j)*SIN(td2)+hh2oav_c2(j)*COS(td2)+ c $ hh2oav_s3(j)*SIN(td3)+hh2oav_c3(j)*COS(td3)+ c $ hh2oav_s4(j)*SIN(td4)+hh2oav_c4(j)*COS(td4) c-------------------------------------------------------------------- ! Water Vapor ! h2o_00=h2oav_00(j)+ $ h2oav_s1(j)*SIN(td1)+ h2oav_c1(j)*COS(td1)+ $ h2oav_s2(j)*SIN(td2)+ h2oav_c2(j)*COS(td2)+ $ h2oav_s3(j)*SIN(td3)+ h2oav_c3(j)*COS(td3)+ $ h2oav_s4(j)*SIN(td4)+ h2oav_c4(j)*COS(td4) ! h2or_m1=h2om1_r00(j)+ $ h2om1_rs1(j)*SIN(td1)+ h2om1_rc1(j)*COS(td1)+ $ h2om1_rs2(j)*SIN(td2)+ h2om1_rc2(j)*COS(td2)+ $ h2om1_rs3(j)*SIN(td3)+ h2om1_rc3(j)*COS(td3)+ $ h2om1_rs4(j)*SIN(td4)+ h2om1_rc4(j)*COS(td4) h2oi_m1=h2om1_i00(j)+ $ h2om1_is1(j)*SIN(td1)+ h2om1_ic1(j)*COS(td1)+ $ h2om1_is2(j)*SIN(td2)+ h2om1_ic2(j)*COS(td2)+ $ h2om1_is3(j)*SIN(td3)+ h2om1_ic3(j)*COS(td3)+ $ h2om1_is4(j)*SIN(td4)+ h2om1_ic4(j)*COS(td4) ! h2or_m2=h2om2_r00(j)+ $ h2om2_rs1(j)*SIN(td1)+ h2om2_rc1(j)*COS(td1)+ $ h2om2_rs2(j)*SIN(td2)+ h2om2_rc2(j)*COS(td2)+ $ h2om2_rs3(j)*SIN(td3)+ h2om2_rc3(j)*COS(td3)+ $ h2om2_rs4(j)*SIN(td4)+ h2om2_rc4(j)*COS(td4) h2oi_m2=h2om2_i00(j)+ $ h2om2_is1(j)*SIN(td1)+ h2om2_ic1(j)*COS(td1)+ $ h2om2_is2(j)*SIN(td2)+ h2om2_ic2(j)*COS(td2)+ $ h2om2_is3(j)*SIN(td3)+ h2om2_ic3(j)*COS(td3)+ $ h2om2_is4(j)*SIN(td4)+ h2om2_ic4(j)*COS(td4) ! h2or_m3=h2om3_r00(j)+ $ h2om3_rs1(j)*SIN(td1)+ h2om3_rc1(j)*COS(td1)+ $ h2om3_rs2(j)*SIN(td2)+ h2om3_rc2(j)*COS(td2)+ $ h2om3_rs3(j)*SIN(td3)+ h2om3_rc3(j)*COS(td3)+ $ h2om3_rs4(j)*SIN(td4)+ h2om3_rc4(j)*COS(td4) h2oi_m3=h2om3_i00(j)+ $ h2om3_is1(j)*SIN(td1)+ h2om3_ic1(j)*COS(td1)+ $ h2om3_is2(j)*SIN(td2)+ h2om3_ic2(j)*COS(td2)+ $ h2om3_is3(j)*SIN(td3)+ h2om3_ic3(j)*COS(td3)+ $ h2om3_is4(j)*SIN(td4)+ h2om3_ic4(j)*COS(td4) ! h2or_m4=h2om4_r00(j)+ $ h2om4_rs1(j)*SIN(td1)+ h2om4_rc1(j)*COS(td1)+ $ h2om4_rs2(j)*SIN(td2)+ h2om4_rc2(j)*COS(td2)+ $ h2om4_rs3(j)*SIN(td3)+ h2om4_rc3(j)*COS(td3)+ $ h2om4_rs4(j)*SIN(td4)+ h2om4_rc4(j)*COS(td4) h2oi_m4=h2om4_i00(j)+ $ h2om4_is1(j)*SIN(td1)+ h2om4_ic1(j)*COS(td1)+ $ h2om4_is2(j)*SIN(td2)+ h2om4_ic2(j)*COS(td2)+ $ h2om4_is3(j)*SIN(td3)+ h2om4_ic3(j)*COS(td3)+ $ h2om4_is4(j)*SIN(td4)+ h2om4_ic4(j)*COS(td4) !-------------------------------------------------------------------- ! scale of Height: wv(1000) and wv(300) ! hh2o_00=hh2oav_00(j)+ $ hh2oav_s1(j)*SIN(td1)+hh2oav_c1(j)*COS(td1)+ $ hh2oav_s2(j)*SIN(td2)+hh2oav_c2(j)*COS(td2)+ $ hh2oav_s3(j)*SIN(td3)+hh2oav_c3(j)*COS(td3)+ $ hh2oav_s4(j)*SIN(td4)+hh2oav_c4(j)*COS(td4) ! hh2or_m1=hh2om1_r00(j)+ $ hh2om1_rs1(j)*SIN(td1)+ hh2om1_rc1(j)*COS(td1)+ $ hh2om1_rs2(j)*SIN(td2)+ hh2om1_rc2(j)*COS(td2)+ $ hh2om1_rs3(j)*SIN(td3)+ hh2om1_rc3(j)*COS(td3)+ $ hh2om1_rs4(j)*SIN(td4)+ hh2om1_rc4(j)*COS(td4) hh2oi_m1=hh2om1_i00(j)+ $ hh2om1_is1(j)*SIN(td1)+ hh2om1_ic1(j)*COS(td1)+ $ hh2om1_is2(j)*SIN(td2)+ hh2om1_ic2(j)*COS(td2)+ $ hh2om1_is3(j)*SIN(td3)+ hh2om1_ic3(j)*COS(td3)+ $ hh2om1_is4(j)*SIN(td4)+ hh2om1_ic4(j)*COS(td4) ! hh2or_m2=hh2om2_r00(j)+ $ hh2om2_rs1(j)*SIN(td1)+ hh2om2_rc1(j)*COS(td1)+ $ hh2om2_rs2(j)*SIN(td2)+ hh2om2_rc2(j)*COS(td2)+ $ hh2om2_rs3(j)*SIN(td3)+ hh2om2_rc3(j)*COS(td3)+ $ hh2om2_rs4(j)*SIN(td4)+ hh2om2_rc4(j)*COS(td4) hh2oi_m2=hh2om2_i00(j)+ $ hh2om2_is1(j)*SIN(td1)+ hh2om2_ic1(j)*COS(td1)+ $ hh2om2_is2(j)*SIN(td2)+ hh2om2_ic2(j)*COS(td2)+ $ hh2om2_is3(j)*SIN(td3)+ hh2om2_ic3(j)*COS(td3)+ $ hh2om2_is4(j)*SIN(td4)+ hh2om2_ic4(j)*COS(td4) ! hh2or_m3=hh2om3_r00(j)+ $ hh2om3_rs1(j)*SIN(td1)+ hh2om3_rc1(j)*COS(td1)+ $ hh2om3_rs2(j)*SIN(td2)+ hh2om3_rc2(j)*COS(td2)+ $ hh2om3_rs3(j)*SIN(td3)+ hh2om3_rc3(j)*COS(td3)+ $ hh2om3_rs4(j)*SIN(td4)+ hh2om3_rc4(j)*COS(td4) hh2oi_m3=hh2om3_i00(j)+ $ hh2om3_is1(j)*SIN(td1)+ hh2om3_ic1(j)*COS(td1)+ $ hh2om3_is2(j)*SIN(td2)+ hh2om3_ic2(j)*COS(td2)+ $ hh2om3_is3(j)*SIN(td3)+ hh2om3_ic3(j)*COS(td3)+ $ hh2om3_is4(j)*SIN(td4)+ hh2om3_ic4(j)*COS(td4) ! hh2or_m4=hh2om4_r00(j)+ $ hh2om4_rs1(j)*SIN(td1)+ hh2om4_rc1(j)*COS(td1)+ $ hh2om4_rs2(j)*SIN(td2)+ hh2om4_rc2(j)*COS(td2)+ $ hh2om4_rs3(j)*SIN(td3)+ hh2om4_rc3(j)*COS(td3)+ $ hh2om4_rs4(j)*SIN(td4)+ hh2om4_rc4(j)*COS(td4) hh2oi_m4=hh2om4_i00(j)+ $ hh2om4_is1(j)*SIN(td1)+ hh2om4_ic1(j)*COS(td1)+ $ hh2om4_is2(j)*SIN(td2)+ hh2om4_ic2(j)*COS(td2)+ $ hh2om4_is3(j)*SIN(td3)+ hh2om4_ic3(j)*COS(td3)+ $ hh2om4_is4(j)*SIN(td4)+ hh2om4_ic4(j)*COS(td4) ! !-------------------------------------------------------------------- ! do i=1,igit dlong1 = 2.*pi*float(i-1)/float(igit) dlong2 = 2.*dlong1 dlong3 = 3.*dlong1 dlong4 = 4.*dlong1 ! ! H2O data ! rh2o(i,j)=h2o_00 $ +h2or_m1*COS(dlong1)+h2oi_m1*SIN(dlong1)+ $ h2or_m2*COS(dlong2)+h2oi_m2*SIN(dlong2)+ $ h2or_m3*COS(dlong3)+h2oi_m3*SIN(dlong3)+ $ h2or_m4*COS(dlong4)+h2oi_m4*SIN(dlong4) H_wv(i,j)=hh2o_00 $ +hh2or_m1*COS(dlong1)+hh2oi_m1*SIN(dlong1)+ $ hh2or_m2*COS(dlong2)+hh2oi_m2*SIN(dlong2)+ $ hh2or_m3*COS(dlong3)+hh2oi_m3*SIN(dlong3)+ $ hh2or_m4*COS(dlong4)+hh2oi_m4*SIN(dlong4) if(rh2o(i,j).le.0.00002)rh2o(i,j)=0.00002 if(H_wv(i,j).ge.3.)H_wv(i,j)=3. ! print*, i,j,rh2o(i,j),H_wv(i,j) ! ! print*, j,hh2o_00,hh2or_m1,hh2oi_m1,hh2or_m2,hh2oi_m2, ! $ hh2or_m3,hh2oi_m3,hh2or_m4,hh2oi_m4 ! c----------------------- UKMO data ip1=i+1 ip2=i+2 im1=i-1 im2=i-2 c--------------------------- if(i.eq.1) then im1=igit im2=igit-1 end if if(i.eq.2)im2=igit c--------------------------- Smoothing the lower bounrady conditions if(i.eq.igit) then ip1=1 ip2=2 end if if(i.eq.igit-1)ip2=1 T_1000(j,i)= (tp_ave(ip2,j)+tp_ave(ip1,j)+tp_ave(i,j)+ $ tp_ave(im1,j)+tp_ave(im2,j))/5. G_1000(j,i)= 9.81*(gh_ave(ip2,j)+gh_ave(ip1,j)+gh_ave(i,j)+ $ gh_ave(im1,j)+gh_ave(im2,j))/5. c T_1000(j,i)= tp_ave(i,j) c G_1000(j,i)= 9.81*gh_ave(i,j) end do 200 continue rb =float(nb) rgit=float(igit) c c.....GEOMETRIC PARAMETERS: c this gives altitude of highest level z_48 =135000. !!!recommended!!! dz = z_48/(48.-0.5) defi = 180./rb ! latitudinal increment dy = a*pi/rb dela = 360./rgit ! longitudinal increment c-------------------- zunten is the lowest level of dynamic fields c----- (!!!vertical wind is calculated at ground, not at zunten!!!) zunten=dz*0.5 dzn = z_48/(float(ngit)-0.5) !!! ngit - dzi in COMMON-BLOCK zunteni=dzn*0.5 C------------------------------- (g/sm^3) dich0 =1.273e-3 rhos =1.273 RgSGS=8.31441e7 RgSI =8314.41 C_p =1005. G =9.81 Rf_c =0.4 co2pmv=360. H =7000. skh=7. c-------------------------------------------------- time step nstep =24*ntime ! steps/day dt =3600./float(ntime) dt2 =dt/2. print 1045,dt 1045 format(10x,'one time step = ',f6.1,'sec',//) do 1 j=1,nb+1 gphi = float(j-1)*defi - 87.5 ! -87.5,...,92.5 North fi =(float(j)-19.5)*defi/180.*pi ! -92.5,...,87.5 North fi2 =(float(j)-19.0)*defi/180.*pi ! -90.0,...,90.0 North c-------------------------------------------------------------------- sinfi(j)=sin(fi) cosfi(j)=cos(fi) cosf2(j)=cos(fi2) sinf2(j)=sin(fi2) tgfia(j)=tan(fi)/a cor(j)=2.*7.292116e-5*sinfi(j) dx(j)=2.*pi/rgit*a*cosfi(j) if(j.eq.nb+1) goto 1 phi(j)=gphi*pi/180. 1 continue cosf2(nb+1)=0. cosf2(1)=0. fi =(float(nb+2)-19.5)*defi/180.*pi ! 92.5 cosfi(nb+2)=cos(fi) c.......................................NEWTONIAN-COOLING-COEFFICIENTS dmax=100. d0=dmax/3. do k=1,kgit z(k)=zunten/1000.+float(k-1)*dz/1000. zw = z(k)*1000.-dz/2. rou(k) = exp(-z(k)*1000./h) row(k) = exp(-zw /h) c........................................... PROFILE OF NEWTONIAN COOLING: alfa(k)=alfaNC(k)*1.e-6 c alfa(k)=(2.3+1.7*tanh((zu-35000.)/8000.)) *1.e-6 ! schoeberl'78 c alfa(k)=(1.5+1.0*tanh((zu-35000.)/7000.)) *1.e-6 ! holton'76 c.........................analytical profile of the eddy diffusion coefficient IF(z(k).GE.90.)geddy0(k)= $ dmax *EXP(-(z(k)-90.)/20.*(z(k)-90.)/20.) IF(z(k).LT.90.)geddy0(k)= $ (dmax-d0)*EXP(-(z(k)-90.)/20.*(z(k)-90.)/20.)+ $ d0 *EXP( (z(k)-90.)/20.) end do c...................................................................ION-DRAG: c if (miond.gt.0) then open (10, FILE='s2_a0_grads.dx' , $ form='unformatted', * access='direct', status='old',recl=nb*23) open (30, FILE='s1u_a0_grads.dx' , $ form='unformatted', * access='direct', status='old',recl=nb*23) open (40, FILE='lh_a0_Mar.dx' , $ form='unformatted', * access='direct', status='old',recl=nb*10) Brayl0=5.E-7 vrAD=PI/180. PI2=2.*PI do k=1,kgit do j=1,nb do i=1,igit drag0(i,j,k,1) = Brayl0 drag0(i,j,k,2) = Brayl0 drag0(i,j,k,3) = 0. end do end do end do c--------------------------------- averaging over time (day) do k=1,23 do j=1,nb s1u_a0(j,k) = 0. s2_a0 (j,k) = 0. end do end do do l=1,24 READ(10,REC=l) s2_a READ(30,REC=l) s1u_a do k=1,23 do j=1,nb s1u_a0(j,k)=s1u_a0(j,k)+s1u_a(j,k)/24. s2_a0 (j,k)=s2_a0 (j,k)+s2_a (j,k)/24. enddo enddo enddo do k=1,10 do j=1,nb ht_a0(j,k) = 0. end do end do c------------------- zonally averaged LH (averaged over January 31 days) do l=1,248 READ(40,REC=l) ht_a do k=1,10 do j=1,nb ht_a0(j,k)=ht_a0(j,k)+ht_a(j,k)/248. enddo enddo enddo c------------------- SPW1 and SPW2 in the LH do k=1,10 do j=1,nb do i=1,igit rlon=float(i-1)*5.625 heat_LH0(i,j,k)=ht_a0(j,k) * +ht_m1(3,j,k,1)*cos(1.*(rlon-ht_m1(3,j,k,3))*VRAD) * +ht_m2(3,j,k,1)*cos(2.*(rlon-ht_m2(3,j,k,3))*VRAD) * +ht_m3(3,j,k,1)*cos(3.*(rlon-ht_m3(3,j,k,3))*VRAD) * +ht_m4(3,j,k,1)*cos(4.*(rlon-ht_m4(3,j,k,3))*VRAD) end do end do end do c-------------------------- Zonally averaged ion drag do k=1,kgit-37 do j=1,nb sinphi=sin(phi(j)) do i =1,igit drag0(i,j,k+37,1)=drag0(i,j,k+37,1)+s1u_a0(j,k) c----------------------------------------- ion drag during the opolar winter IF (drag0(i,j,k+37,1).LT.5.e-7) drag0(i,j,k+37,1)=5.E-7 IF (drag0(i,j,k+37,1).LT.5.e-6.and.k+37.gt.44) & drag0(i,j,k+37,1)=5.E-6 drag0(i,j,k+37,2)=drag0(i,j,k+37,2)+s1u_a0(j,k)* & 4.*sinphi*sinphi/(1.+3.*sinphi*sinphi) drag0(i,j,k+37,3)=s2_a0(j,k) drag(i,j,k+37,1)=drag0(i,j,k+37,1) drag(i,j,k+37,2)=drag0(i,j,k+37,2) drag(i,j,k+37,3)=drag0(i,j,k+37,3) enddo enddo end do c end if c................................GLOBAL MEAN TEMPERATURE tz(kgit) c----------------- initial setting, will be improved in strobel.f do k=1,kgit tz(k)=0. do j=1,nb do i=1,igit tz(k)=tz(k)+tn1(j,k,i)/rgit/rb enddo enddo end do C................................... DENSITY and related PROFILES c----------------- initial setting, will be improved in strobel.f c---------- T(0), pressure at the lower boundary is 1.013e6 - SGS t0=1.013e6*rm(1)/RgSGS/dich0 do k=1,kgit Hsclkm=RgSGS/rm(k)*tz(k)*1.e-7/9.81 dichte(k)=dich0*rou(k)*rm(k)/rm(1)*t0/tz(k) c------------------------------------- xmt will be used in grav.f xmt(k)=rm(k)/rm(1)*t0/tz(k) c------------------------------------------ for UV radiation: c c----------------- Cp(z) - approx. 1.e7 (SGS) in lower atmosphere c cp_z=RgSGS/rm(k)*(7./2.*(1.-roxvmr(k))+5./2.*roxvmr(k)) dichcp(k)=dichte(k)*cp_z sunris(k)=acos(6373./(6373.+z(k)))*180./pi + 90. do j=1,nb c c------------- initial setting, will be recalculated in strobel.f c do i=1,igit xnumo3(j,k,i)=1.013 *rou(k)*ozppmv(j,k,i)/1.38066e-16/tz(k) end do xnumn2(j,k)=1.013e6*rou(k)*rn2vmr( k)/1.38066e-16/tz(k) xnumo2(j,k)=1.013e6*rou(k)*ro2vmr( k)/1.38066e-16/tz(k) xnumox(j,k)=1.013e6*rou(k)*roxvmr( k)/1.38066e-16/tz(k) sox(j,k)=Hsclkm*rm(k)/16.*xnumox(j,k)*1.e5 so2(j,k)=Hsclkm*rm(k)/32.*xnumo2(j,k)*1.e5 sn2(j,k)=Hsclkm*rm(k)/28.*xnumn2(j,k)*1.e5 end do end do c---- M/M0*T0/T xmt(kgit+1)=(8.*xmt(kgit)-3.*xmt(kgit-1))/5. c c....................................................H2O: c ratio=28.9644/18.016 !mol_luft/mol_h2o c------------------------ this is the default value at ground c--------- new value Forbes and Garrett (1979) c ppm=0.0175 c ppm=0.02 !max.massemmisch.verh.(=0.032ppv) c ppm=0.03 !max.massemmisch.verh.(=0.048ppv) c c c2_h2o=2.e-6 do 12 i=1,igit do 12 j=1,nb D_1000=1.e6*rm(1)/RgSGS/T_1000(j,i) c D_1000=1.e6*rm(1)/RgSGS/273. c gphi = -87.5+float(j-1)*5. ! -87.5,...,87.5 North c---------------------ANALITICAL H2O(phi) c c1_h2o(i,j)= ppm*exp (-gphi*gphi/1600.) c a1h2o(i,j) = exp (-gphi*gphi/1600.)*ppm*dich0 *1.e5 c a1h2o(i,j) = exp (-gphi*gphi/1600.)*ppm*D_1000*1.e5 c------------------- rh2o(j) is the empirical (NCEP/NCAR) distribution c D_1000 (g/cm^3), 1.e5 (km - cm), a1h2o (g/cm^2/km) c------------------NCEP/NCAR H2O DATA--------------------------------- ! ! NCEP/NCAR H2O 3D DATA (NEW!) ! a1h2o(i,j)= rh2o(i,j)*D_1000*1.e5 ! a1h2o(i,j)= rh2o(j)*D_1000*1.e5 12 continue ! !------------------- factors for h2o column !------------------- will be used in ircool_newCO2.f and strobel.f do 13 kk=1,4 skhnir(kk) = skh/(1.+xncir(kk)) do i=1,igit do j=1,nb e1nir (kk,i,j) = -1./skhnir(kk) - 1./H_wv(i,j) end do end do 13 continue do 14 kk=1,9 skhnuv(kk) = skh/(1.+dkli(kk)) do i=1,igit do j=1,nb expnuv(kk,i,j) = -1./skhnuv(kk) - 1./H_wv(i,j) end do end do 14 continue c---------------------------------- new H2O heating??????? c do 15 kk=1,53 c H_h2o_2(kk) = skh/(1.+cm_h2o(kk)) c do j=1,nb c H_h2o_1(kk,j) = 1./(1./H_h2o_2(kk)+1./H_wv(j)) c end do c----------------- scaling factors for each interval c s1_h2o(kk)=(1013./p0_h2o(kk))**cm_h2o(kk) c 15 continue return end