wrf-fire /wrfv2_fire/phys/module_bl_boulac.F

Language Fortran 77 Lines 922
MD5 Hash 0590a77da7a2c9d54b43c8e91ec7c0bb Estimated Cost $15,608 (why?)
Repository git://github.com/jbeezley/wrf-fire.git View Raw File View Project SPDX
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
MODULE module_bl_boulac

!USE module_model_constants
    

!------------------------------------------------------------------------
!    Calculation of the tendency due to momentum, heat 
!    and moisture turbulent fluxes follwing the approach 
!    of Bougeault and Lacarrere, 1989 (MWR, 117, 1872-1890).
!    The scheme computes a prognostic ecuation for TKE and derives 
!    dissipation and turbulent coefficients using length scales.
!    
!    Subroutine written by Alberto Martilli, CIEMAT, Spain,
!    e-mail:alberto_martilli@ciemat.es
!    August 2006.
!------------------------------------------------------------------------
! IN THIS VERSION TKE IS NOT ADVECTED!!!!
! TO BE CHANGED IN THE FUTURE
! 
! -----------------------------------------------------------------------
!  Constant used in the module
!  ck_b=constant used in the compuation of diffusion coefficients
!  ceps_b=constant used inthe computation of dissipation
!  temin= minimum value allowed for TKE
!  vk=von karman constant
! -----------------------------------------------------------------------

      real ck_b,ceps_b,vk,temin    ! constant for Bougeault and Lacarrere    
      parameter(ceps_b=1/1.4,ck_b=0.4,temin=0.0001,vk=0.4) ! impose minimum values for tke similar to those of MYJ
! -----------------------------------------------------------------------     


   CONTAINS
 
      subroutine boulac(frc_urb2d,idiff,flag_bep,dz8w,dt,u_phy,v_phy   & 
                      ,th_phy,rho,qv_curr,hfx                                  &
                      ,qfx,ustar,cp,g                                          &
                      ,rublten,rvblten,rthblten                                &
                      ,rqvblten,rqcblten                        & 
                      ,tke,dlk,wu,wv,wt,wq,exch_h,exch_m,pblh        &
                      ,a_u_bep,a_v_bep,a_t_bep,a_q_bep          &
                      ,a_e_bep,b_u_bep,b_v_bep                  &
                      ,b_t_bep,b_q_bep,b_e_bep,dlg_bep          &
                      ,dl_u_bep,sf_bep,vl_bep                   &                 
                      ,ids,ide, jds,jde, kds,kde                &
                      ,ims,ime, jms,jme, kms,kme                &
                      ,its,ite, jts,jte, kts,kte)                    

      implicit none



!-----------------------------------------------------------------------
!     Input
!------------------------------------------------------------------------
   INTEGER::                        ids,ide, jds,jde, kds,kde,  &
                                    ims,ime, jms,jme, kms,kme,  &
                                    its,ite, jts,jte, kts,kte
 
   integer, INTENT(IN) :: idiff
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   DZ8W      !vertical resolution       
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   qv_curr   !moisture  
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   RHO       !air density
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   TH_PHY    !potential temperature
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   U_PHY     !x-component of wind
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::   V_PHY     !y-component of wind
   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::   ustar              !friction velocity
   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::   hfx                !sensible heat flux (W/m2) at surface 
   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::   qfx                !moisture flux at surface
   real,  INTENT(IN   )    :: g,cp                                              !gravity and Cp
   REAL, INTENT(IN )::   DT                                                      ! Time step

   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::   FRC_URB2D          !fraction cover urban
   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)    ::   PBLH          !PBL height
!
! variable added for urban
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::a_u_bep        ! Implicit component for the momemtum in X-direction
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::a_v_bep        ! Implicit component for the momemtum in Y-direction
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::a_t_bep        ! Implicit component for the Pot. Temp.
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::a_q_bep        ! Implicit component for Moisture
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::a_e_bep        ! Implicit component for the TKE
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::b_u_bep        ! Explicit component for the momemtum in X-direction
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::b_v_bep        ! Explicit component for the momemtum in Y-direction
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::b_t_bep        ! Explicit component for the Pot. Temp.
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::b_q_bep        ! Explicit component for Moisture
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::b_e_bep        ! Explicit component for the TKE

   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT)    ::dlg_bep        ! Height above ground (L_ground in formula (24) of the BLM paper). 
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::dl_u_bep        ! Length scale (lb in formula (22) ofthe BLM paper).
! urban surface and volumes        
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::sf_bep           ! surface of the urban grid cells
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN   )    ::vl_bep             ! volume of the urban grid cells
   LOGICAL, INTENT(IN) :: flag_bep                                             !flag for BEP
                                                           
!
!-----------------------------------------------------------------------
!     Local, carried on from one timestep to the other
!------------------------------------------------------------------------
!      real, save, allocatable, dimension (:,:,:)::TKE  ! Turbulent kinetic energy
      real,  dimension (ims:ime, kms:kme, jms:jme)  ::th_0 ! reference state for potential temperature

!------------------------------------------------------------------------
!     Output
!------------------------------------------------------------------------   
        real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  exch_h ! exchange coefficient for heat
        real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  exch_m ! exchange coefficient for momentum
        real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(INOUT   )    ::  tke  ! Turbulence Kinetic Energy 
        real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  wu  ! Turbulent flux of momentum (x) 
        real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  wv  ! Turbulent flux of momentum (y) 
        real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  wt  ! Turbulent flux of temperature
        real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  wq  ! Turbulent flux of water vapor
        real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::  dlk  ! Turbulent flux of water vapor
! only if idiff not equal 1:
        REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::   RUBLTEN  !tendency for U_phy
        REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::   RVBLTEN  !tendency for V_phy
        REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::   RTHBLTEN !tendency for TH_phy
        REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::   RQVBLTEN !tendency for QV_curr
        REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT   )    ::   RQCBLTEN !tendency for QV_curr

!--------------------------------------------------------------
! Local
!--------------------------------------------------------------
! 1D array used for the input and output of the routine boulac1D

      real z1D(kms:kme)               ! vertical coordinates (faces of the grid)
      real dz1D(kms:kme)              ! vertical resolution
      real u1D(kms:kme)                 ! wind speed in the x directions
      real v1D(kms:kme)                 ! wind speed in the y directions
      real th1D(kms:kme)                ! potential temperature
      real q1D(kms:kme)                 ! potential temperature
      real rho1D(kms:kme)               ! air density
      real rhoz1D(kms:kme)            ! air density at the faces
      real tke1D(kms:kme)               ! air pressure
      real th01D(kms:kme)               ! reference potential temperature
      real dlk1D(kms:kme)               ! dlk
      real dls1D(kms:kme)               ! dls
      real exch1D(kms:kme)            ! exch
      real sf1D(kms:kme)              ! surface of the grid cells
      real vl1D(kms:kme)                ! volume of the  grid cells
      real a_u1D(kms:kme)               ! Implicit component of the momentum sources or sinks in the X-direction
      real a_v1D(kms:kme)               ! Implicit component of the momentum sources or sinks in the Y-direction
      real a_t1D(kms:kme)               ! Implicit component of the heat sources or sinks
      real a_q1D(kms:kme)               ! Implicit component of the moisture sources or sinks
      real a_e1D(kms:kme)               ! Implicit component of the TKE sources or sinks
      real b_u1D(kms:kme)               ! Explicit component of the momentum sources or sinks in the X-direction
      real b_v1D(kms:kme)               ! Explicit component of the momentum sources or sinks in the Y-direction
      real b_t1D(kms:kme)               ! Explicit component of the heat sources or sinks
      real b_q1D(kms:kme)               ! Explicit component of the moisture sources or sinks
      real b_e1D(kms:kme)               ! Explicit component of the TKE sources or sinks
      real dlg1D(kms:kme)               ! Height above ground (L_ground in formula (24) of the BLM paper). 
      real dl_u1D(kms:kme)              ! Length scale (lb in formula (22) ofthe BLM paper)
      real sh1D(kms:kme)              ! shear
      real bu1D(kms:kme)              ! buoyancy
      real wu1D(kms:kme)              ! turbulent flux of momentum (x component)
      real wv1D(kms:kme)              ! turbulent flux of momentum (y component)
      real wt1D(kms:kme)              ! turbulent flux of temperature
      real wq1D(kms:kme)              ! turbulent flux of water vapor
! local added only for diagnostic output
      real a_e(ims:ime,kms:kme,jms:jme) ! implicit term in TKE
      real b_e(ims:ime,kms:kme,jms:jme) ! explicit term in TKE
      real bu(ims:ime,kms:kme,jms:jme) ! buoyancy term in TKE
      real sh(ims:ime,kms:kme,jms:jme) ! shear term in TKE
      real wrk(ims:ime) ! working array
      integer ix,iy,iz,id,iz_u,iw_u,ig,ir_u,ix1,iy1
      real ufrac_int                                              ! urban fraction     
      real vect,time_tke,hour,zzz
      real ustarf
      real summ1,summ2,summ3
      save time_tke,hour
!
! store reference state for potential temperature
! and initialize tke  
!    

!here I fix the value of the reference state equal to the value of the potnetial temperature
! the only use of this variable in the code is to compute the paramter BETA = g/th0
! I fix it to 300K. 
      
        do ix=its,ite
        do iy=jts,jte        
        do iz=kts,kte
!         th_0(ix,iz,iy)=th_phy(ix,iz,iy)
         th_0(ix,iz,iy)=300.
        enddo
        enddo
        enddo
        
         bu1d=0.
         sh1d=0.         
         b_e1d=0.
         b_u1d=0.
         b_v1d=0.
         b_t1d=0.
         b_q1d=0.
         a_e1d=0.
         a_u1d=0.
         a_v1d=0.
         a_t1d=0.
         a_q1d=0.
       z1D=0.               ! vertical coordinates (faces of the grid)
       dz1D=0.              ! vertical resolution
       u1D =0.                ! wind speed in the x directions
       v1D =0.                ! wind speed in the y directions
       th1D=0.                ! potential temperature
       q1D=0.                 ! potential temperature
       rho1D=0.               ! air density
       rhoz1D=0.            ! air density at the faces
       tke1D =0.              ! air pressure
       th01D =0.              ! reference potential temperature
       dlk1D =0.              ! dlk
       dls1D =0.              ! dls
       exch1D=0.            ! exch
       sf1D  =1.            ! surface of the grid cells
       vl1D  =1.             ! volume of the  grid cells
       a_u1D =0.              ! Implicit component of the momentum sources or sinks in the X-direction
       a_v1D =0.              ! Implicit component of the momentum sources or sinks in the Y-direction
       a_t1D =0.              ! Implicit component of the heat sources or sinks
       a_q1D =0.              ! Implicit component of the moisture sources or sinks
       a_e1D =0.              ! Implicit component of the TKE sources or sinks
       b_u1D =0.              ! Explicit component of the momentum sources or sinks in the X-direction
       b_v1D =0.              ! Explicit component of the momentum sources or sinks in the Y-direction
       b_t1D =0.              ! Explicit component of the heat sources or sinks
       b_q1D =0.              ! Explicit component of the moisture sources or sinks
       b_e1D =0.              ! Explicit component of the TKE sources or sinks
       dlg1D =0.              ! Height above ground (L_ground in formula (24) of the BLM paper). 
       dl_u1D=0.              ! Length scale (lb in formula (22) ofthe BLM paper)
       sh1D  =0.            ! shear
       bu1D  =0.            ! buoyancy
       wu1D  =0.            ! turbulent flux of momentum (x component)
       wv1D  =0.            ! turbulent flux of momentum (y component)
       wt1D =0.              ! turbulent flux of temperature
       wq1D =0.             ! turbulent flux of water vapor
! local added only for diagnostic output
         
! loop over the columns. 
! put variables in 1D temporary arrays
!     

       do ix=its,ite
       do iy=jts,jte
         z1d(kts)=0.
         do iz= kts,kte
	  u1D(iz)=u_phy(ix,iz,iy)
	  v1D(iz)=v_phy(ix,iz,iy)
	  th1D(iz)=th_phy(ix,iz,iy)
          q1D(iz)=qv_curr(ix,iz,iy)
          tke1D(iz)=tke(ix,iz,iy)
	  rho1D(iz)=rho(ix,iz,iy)	   
	  th01D(iz)=th_0(ix,iz,iy)	  
          dz1D(iz)=dz8w(ix,iz,iy)
          z1D(iz+1)=z1D(iz)+dz1D(iz)
         enddo
        rhoz1D(kts)=rho1D(kts)
        do iz=kts+1,kte
         rhoz1D(iz)=(rho1D(iz)*dz1D(iz-1)+rho1D(iz-1)*dz1D(iz))/(dz1D(iz-1)+dz1D(iz))
        enddo
        rhoz1D(kte+1)=rho1D(kte)
        if(flag_bep)then
         do iz=kts,kte          
          a_e1D(iz)=a_e_bep(ix,iz,iy)
          b_e1D(iz)=b_e_bep(ix,iz,iy)
          dlg1D(iz)=(z1D(iz)+z1D(iz+1))/2.*(1.-frc_urb2d(ix,iy))+dlg_bep(ix,iz,iy)*frc_urb2d(ix,iy)
          dl_u1D(iz)=dl_u_bep(ix,iz,iy)
          if((1.-frc_urb2d(ix,iy)).lt.1.)dl_u1D(iz)=dl_u1D(iz)/frc_urb2d(ix,iy)
          vl1D(iz)=vl_bep(ix,iz,iy)
          sf1D(iz)=sf_bep(ix,iz,iy)
         enddo
         ufrac_int=frc_urb2d(ix,iy)
         sf1D(kte+1)=sf_bep(ix,1,iy)
        else
         do iz=kts,kte          
          a_e1D(iz)=0.        
          b_e1D(iz)=0.
          dlg1D(iz)=(z1D(iz)+z1D(iz+1))/2.
          dl_u1D(iz)=0.
          vl1D(iz)=1.
          sf1D(iz)=1.
         enddo
         ufrac_int=0.
         sf1D(kte+1)=1.
        endif
! call the routine that will solve the turbulence in 1D for tke

         call boulac1D(ix,iy,ufrac_int,kms,kme,kts,kte,dz1d,z1D,dt,u1D,v1D,th1D,rho1D,rhoz1D,q1D,th01D,&
                   tke1D,ustar(ix,iy),hfx(ix,iy),qfx(ix,iy),cp,g,        & 
                   a_e1D,b_e1D,                              & 
                   dlg1D,dl_u1D,sf1D,vl1D,dlk1D,dls1D,exch1D,sh1D,bu1D)                    

         call pbl_height(kms,kme,kts,kte,dz1d,z1d,th1D,q1D,pblh(ix,iy))


! store turbulent exchange coefficients, TKE, and other variables
         
         do iz= kts,kte
          a_e(ix,iz,iy)=a_e1D(iz)
          b_e(ix,iz,iy)=b_e1D(iz)
          if(flag_bep)then
          dlg_bep(ix,iz,iy)=dlg1D(iz)
          endif
          tke(ix,iz,iy)=tke1D(iz)
          dlk(ix,iz,iy)=dlk1D(iz)
          sh(ix,iz,iy)=sh1D(iz)
          bu(ix,iz,iy)=bu1D(iz)
          exch_h(ix,iz,iy)=exch1D(iz)
          exch_m(ix,iz,iy)=exch1D(iz)
         enddo

         if(idiff.ne.1)then

! estimate the tendencies

        if(flag_bep)then         
         do iz=kts,kte          
          a_t1D(iz)=a_t_bep(ix,iz,iy)
          b_t1D(iz)=b_t_bep(ix,iz,iy)
          a_u1D(iz)=a_u_bep(ix,iz,iy)
          b_u1D(iz)=b_u_bep(ix,iz,iy)
          a_v1D(iz)=a_v_bep(ix,iz,iy)
          b_v1D(iz)=b_v_bep(ix,iz,iy)
          a_q1D(iz)=a_q_bep(ix,iz,iy)
          b_q1D(iz)=b_q_bep(ix,iz,iy)
         enddo
        else
         do iz=kts,kte          
          a_t1D(iz)=0.         
          b_t1D(iz)=0.
          a_u1D(iz)=0.        
          b_u1D(iz)=0.
          a_v1D(iz)=0.         
          b_v1D(iz)=0.
          a_q1D(iz)=0.        
          b_q1D(iz)=0.
         enddo
          b_t1D(1)=hfx(ix,iy)/dz1D(1)/rho1D(1)/cp         
          b_q1D(1)=qfx(ix,iy)/dz1D(1)/rho1D(1)
          a_u1D(1)=(-ustar(ix,iy)*ustar(ix,iy)/dz1D(1)/((u1D(1)**2.+v1D(1)**2.)**.5))
          a_v1D(1)=(-ustar(ix,iy)*ustar(ix,iy)/dz1D(1)/((u1D(1)**2.+v1D(1)**2.)**.5))
        endif

! solve diffusion equation for momentum x component          
       call diff(kms,kme,kts,kte,1,1,dt,u1D,rho1D,rhoz1D,exch1D,a_u1D,b_u1D,sf1D,vl1D,dz1D,wu1D)
        
! solve diffusion equation for momentum y component          
       call diff(kms,kme,kts,kte,1,1,dt,v1D,rho1D,rhoz1D,exch1D,a_v1D,b_v1D,sf1D,vl1D,dz1D,wv1D)

! solve diffusion equation for potential temperature        
       call diff(kms,kme,kts,kte,1,1,dt,th1D,rho1D,rhoz1D,exch1D,a_t1D,b_t1D,sf1D,vl1D,dz1D,wt1D)

! solve diffusion equation for water vapor mixing ratio     
       call diff(kms,kme,kts,kte,1,1,dt,q1D,rho1D,rhoz1D,exch1D,a_q1D,b_q1D,sf1D,vl1D,dz1D,wq1D)

! compute the tendencies
                             
         do iz= kts,kte
          rthblten(ix,iz,iy)=rthblten(ix,iz,iy)+(th1D(iz)-th_phy(ix,iz,iy))/dt
          rqvblten(ix,iz,iy)=rqvblten(ix,iz,iy)+(q1D(iz)-qv_curr(ix,iz,iy))/dt
          rublten(ix,iz,iy)=rublten(ix,iz,iy)+(u1D(iz)-u_phy(ix,iz,iy))/dt
          rvblten(ix,iz,iy)=rvblten(ix,iz,iy)+(v1D(iz)-v_phy(ix,iz,iy))/dt  
          wu(ix,iz,iy)=wu1D(iz)
          wv(ix,iz,iy)=wv1D(iz) 
          wt(ix,iz,iy)=wt1D(iz)
          wq(ix,iz,iy)=wq1D(iz)      
        enddo
      endif
   
      enddo  ! iy
      enddo  ! ix

         
         time_tke=time_tke+dt

  
      return
      end subroutine boulac
            

! ===6=8===============================================================72

! ===6=8===============================================================72

      subroutine boulac1D(ix,iy,ufrac_int,kms,kme,kts,kte,dz,z,dt,u,v,th,rho,rhoz,qa,th0,te,    &
                   ustar,hfx,qfx,cp,g,                               & 
                   a_e,b_e,                        & 
                   dlg,dl_u,sf,vl,dlk,dls,exch,sh,bu)                           

! ----------------------------------------------------------------------
! 1D resolution of TKE following Bougeault and Lacarrere
! ----------------------------------------------------------------------

      implicit none

      integer iz,ix,iy

! ----------------------------------------------------------------------
! INPUT:
! ----------------------------------------------------------------------


      integer kms,kme,kts,kte                 
      real z(kms:kme)               ! Altitude above the ground of the cell interfaces.
      real dz(kms:kme)                ! vertical resolution
      real u(kms:kme)                ! Wind speed in the x direction
      real v(kms:kme)                ! Wind speed in the y direction
      real th(kms:kme)                ! Potential temperature
      real rho(kms:kme)                ! Air density
      real g                     ! gravity
      real cp                    !  
      real te(kms:kme)                ! turbulent kinetic energy
      real qa(kms:kme)                ! air humidity
      real th0(kms:kme)               ! Reference potential temperature 
      real dt                    ! Time step
      real ustar                 ! ustar
      real hfx                   ! sensbile heat flux
      real qfx                   ! kinematic latent heat flux
      real sf(kms:kme)              ! surface of the urban grid cells
      real vl(kms:kme)                ! volume of the urban grid cells
      real a_e(kms:kme)               ! Implicit component of the TKE sources or sinks
      real b_e(kms:kme)               ! Explicit component of the TKE sources or sinks
      real dlg(kms:kme)               ! Height above ground (L_ground in formula (24) of the BLM paper). 
      real dl_u(kms:kme)              ! Length scale (lb in formula (22) ofthe BLM paper)
      real ufrac_int             ! urban fraction
! local variables not needed in principle, but that can be used to estimate the budget and turbulent fluxes
 
      real we(kms:kme),dwe(kms:kme)
     
! local variables
      real sh(kms:kme)    ! shear term in TKE eqn.
      real bu(kms:kme)    ! buoyancy term in TKE eqn.
      real td(kms:kme)    ! dissipation term in TKE eqn.
      real exch(kms:kme) ! turbulent diffusion coefficients (defined at the faces)
      real dls(kms:kme)   ! dissipation length scale
      real dlk(kms:kme)   ! length scale used to estimate exch
      real dlu(kms:kme)   ! l_up
      real dld(kms:kme)   ! l_down
      real rhoz(kms:kme) !air density at the faces of the cell                
      real tstar     ! derived from hfx and ustar
      real beta
      real summ1,summ2,summ3,summ4
! interpolate air density at the faces

       

! estimation of tstar

        tstar=-hfx/rho(1)/cp/ustar                                                
         
! first compute values of dlu and dld (length scales up and down). 
       
       call dissip_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,te,dlu,dld,th,th0)
          
!then average them to obtain dls and dlk (length scales for dissipation and eddy coefficients)

       call length_bougeault(ix,iy,kms,kme,kts,kte,dld,dlu,dlg,dl_u,dls,dlk)
            
! compute the turbulent diffusion coefficients exch

       call cdtur_bougeault(ix,iy,kms,kme,kts,kte,te,z,dz,exch,dlk)

! compute source and sink terms in the TKE equation (shear, buoyancy and dissipation)
       
       call tke_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,vl,u,v,th,te,th0,ustar,tstar,exch,dls,td,sh,bu,b_e,a_e,sf,ufrac_int)
        
! solve for tke 
      
       call diff(kms,kme,kts,kte,1,1,dt,te,rho,rhoz,exch,a_e,b_e,sf,vl,dz,we)
     
! avoid negative values for tke
  
       do iz=kts,kte
        if(te(iz).lt.temin) te(iz)=temin
       enddo 
      
       return
       end subroutine boulac1d
! 
! ===6=8===============================================================72

! ===6=8===============================================================72
         subroutine dissip_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,te,dlu,dld,th,th0)
! compute the length scales up and down
         implicit none
         integer kms,kme,kts,kte,iz,izz,ix,iy
         real dzt,zup,beta,zup_inf,bbb,tl,zdo,zdo_sup,zzz,g
         real te(kms:kme),dlu(kms:kme),dld(kms:kme),dz(kms:kme)
         real z(kms:kme),th(kms:kme),th0(kms:kme)

         do iz=kts,kte
          zup=0.
          dlu(iz)=z(kte+1)-z(iz)-dz(iz)/2.
          zzz=0.
          zup_inf=0.
          beta=g/th0(iz)      !Buoyancy coefficient
          do izz=iz,kte-1
           dzt=(dz(izz+1)+dz(izz))/2.
           zup=zup-beta*th(iz)*dzt
           zup=zup+beta*(th(izz+1)+th(izz))*dzt/2.
           zzz=zzz+dzt
           if(te(iz).lt.zup.and.te(iz).ge.zup_inf)then
            bbb=(th(izz+1)-th(izz))/dzt
            if(bbb.ne.0)then
             tl=(-beta*(th(izz)-th(iz))+sqrt( max(0.,(beta*(th(izz)-th(iz)))**2.+2.*bbb*beta*(te(iz)-zup_inf))))/bbb/beta
            else
             if(th(izz).ne.th(iz))then
              tl=(te(iz)-zup_inf)/(beta*(th(izz)-th(iz)))
             else
              tl=0.
             endif
            endif            
            dlu(iz)=zzz-dzt+tl
           endif
           zup_inf=zup
          enddo
                  
          zdo=0.
          zdo_sup=0.
          dld(iz)=z(iz)+dz(iz)/2.
          zzz=0.
          do izz=iz,kts+1,-1
           dzt=(dz(izz-1)+dz(izz))/2.
           zdo=zdo+beta*th(iz)*dzt
           zdo=zdo-beta*(th(izz-1)+th(izz))*dzt/2.
           zzz=zzz+dzt
           if(te(iz).lt.zdo.and.te(iz).ge.zdo_sup)then
            bbb=(th(izz)-th(izz-1))/dzt
            if(bbb.ne.0.)then
             tl=(beta*(th(izz)-th(iz))+sqrt( max(0.,(beta*(th(izz)-th(iz)))**2.+2.*bbb*beta*(te(iz)-zdo_sup))))/bbb/beta
            else
             if(th(izz).ne.th(iz))then
              tl=(te(iz)-zdo_sup)/(beta*(th(izz)-th(iz)))
             else
              tl=0.
             endif
            endif
            
            dld(iz)=zzz-dzt+tl
           endif
           zdo_sup=zdo
          enddo
         enddo
            
                   
         end subroutine dissip_bougeault
!
! ===6=8===============================================================72 
! ===6=8===============================================================72
         subroutine length_bougeault(ix,iy,kms,kme,kts,kte,dld,dlu,dlg,dl_u,dls,dlk)
! compute the length scales for dissipation and turbulent coefficients
         implicit none
         integer kms,kme,kts,kte,iz,ix,iy
         real dlu(kms:kme),dld(kms:kme),dl_u(kms:kme)
         real dls(kms:kme),dlk(kms:kme),dlg(kms:kme)
         
         do iz=kts,kte
          dld(iz)=min(dld(iz),dlg(iz))
          dls(iz)=sqrt(dlu(iz)*dld(iz))
          dlk(iz)=min(dlu(iz),dld(iz))

         if(dl_u(iz).gt.0.)then               
           dls(iz)=1./(1./dls(iz)+1./dl_u(iz))
           dlk(iz)=1./(1./dlk(iz)+1./dl_u(iz))             
          endif
         enddo 
                   
         return
         end subroutine length_bougeault
!

! ===6=8===============================================================72 
! ===6=8===============================================================72

       subroutine cdtur_bougeault(ix,iy,kms,kme,kts,kte,te,z,dz,exch,dlk)
! compute turbulent coefficients
       implicit none
       integer iz,kms,kme,kts,kte,ix,iy
       real te_m,dlk_m
       real te(kms:kme),exch(kms:kme)
       real dz(kms:kme),z(kms:kme)
       real dlk(kms:kme)
       real fact

       exch(kts)=0.

!       do iz=2,nz-1
       do iz=kts+1,kte
        te_m=(te(iz-1)*dz(iz)+te(iz)*dz(iz-1))/(dz(iz)+dz(iz-1))
        dlk_m=(dlk(iz-1)*dz(iz)+dlk(iz)*dz(iz-1))/(dz(iz)+dz(iz-1))
        exch(iz)=ck_b*dlk_m*sqrt(te_m)
!        exch(iz)=max(exch(iz),0.0001)    
        exch(iz)=max(exch(iz),0.1) 
       enddo

       exch(kte+1)=0.1

       return
       end subroutine cdtur_bougeault


! ===6=8===============================================================72
! ===6=8===============================================================72

       subroutine diff(kms,kme,kts,kte,iz1,izf,dt,co,rho,rhoz,cd,aa,bb,sf,vl,dz,fc)


!------------------------------------------------------------------------
!           Calculation of the diffusion in 1D        
!------------------------------------------------------------------------
!  - Input:
!       nz    : number of points
!       iz1   : first calculated point
!       co    : concentration of the variable of interest
!       dz    : vertical levels
!       cd    : diffusion coefficients
!       dtext : external time step
!       rho    : density of the air at the center
!       rhoz   : density of the air at the face
!       itest : if itest eq 1 then update co, else store in a flux array
!  - Output:
!       co :concentration of the variable of interest

!  - Internal:
!       cddz  : constant terms in the equations 
!       dt    : diffusion time step
!       nt    : number of the diffusion time steps
!       cstab : ratio of the stability condition for the time step
!---------------------------------------------------------------------

        implicit none
        integer iz,iz1,izf
        integer kms,kme,kts,kte
        real dt,dzv               
        real co(kms:kme),cd(kms:kme),dz(kms:kme)
        real rho(kms:kme),rhoz(kms:kme)
        real cddz(kms:kme+1),fc(kms:kme),df(kms:kme)
        real a(kms:kme,3),c(kms:kme)
        real sf(kms:kme),vl(kms:kme)
        real aa(kms:kme),bb(kms:kme)
        

! Compute cddz=2*cd/dz  
        
        cddz(kts)=sf(kts)*rhoz(kts)*cd(kts)/dz(kts)
        do iz=kts+1,kte
         cddz(iz)=2.*sf(iz)*rhoz(iz)*cd(iz)/(dz(iz)+dz(iz-1))
        enddo
        cddz(kte+1)=sf(kte+1)*rhoz(kte+1)*cd(kte+1)/dz(kte)

         do iz=kts,iz1-1
          a(iz,1)=0.
          a(iz,2)=1.
          a(iz,3)=0.
          c(iz)=co(iz)
         enddo
          
          do iz=iz1,kte-izf
           dzv=vl(iz)*dz(iz)
           a(iz,1)=-cddz(iz)*dt/dzv/rho(iz)
           a(iz,2)=1+dt*(cddz(iz)+cddz(iz+1))/dzv/rho(iz)-aa(iz)*dt
           a(iz,3)=-cddz(iz+1)*dt/dzv/rho(iz)
           c(iz)=co(iz)+bb(iz)*dt                     
          enddo
          
          do iz=kte-(izf-1),kte
           a(iz,1)=0.
           a(iz,2)=1
           a(iz,3)=0.
           c(iz)=co(iz)
          enddo
           
          call invert (kms,kme,kts,kte,a,c,co)
         
          do iz=kts,iz1 
           fc(iz)=0.
          enddo
                       
          do iz=iz1+1,kte 
           fc(iz)=-(cddz(iz)*(co(iz)-co(iz-1)))/rho(iz)
          enddo
        
!          do iz=1,iz1
!           df(iz)=0.
!          enddo
!          
!          do iz=iz1+1,nz-izf
!           dzv=vl(iz)*dz(iz)
!           df(iz)=+(co(iz-1)*cddz(iz)-co(iz)*(cddz(iz)+cddz(iz+1))+co(iz+1)*cddz(iz+1))/dzv/rho(iz)
!          enddo
!          
!          do iz=nz-izf,nz 
!           df(iz)=0.
!          enddo
                                        
       return
       end subroutine diff

! ===6=8===============================================================72
! ===6=8===============================================================72

       subroutine buoy(ix,iy,g,kms,kme,kts,kte,th,th0,exch,dz,bu,ustar,tstar,ufrac_int)
! compute buoyancy term
       implicit none
       integer kms,kme,kts,kte,iz,ix,iy
       real dtdz1,dtdz2,cdm,dtmdz,g      
       real th(kms:kme),exch(kms:kme),dz(kms:kme),bu(kms:kme)
       real th0(kms:kme),ustar,tstar,ufrac_int
        
!       bu(1)=-ustar*tstar*g/th0(1)*(1.-ufrac_int) 
      bu(kts)=0. 
       
        
       do iz=kts+1,kte-1       
        dtdz1=2.*(th(iz)-th(iz-1))/(dz(iz-1)+dz(iz))
        dtdz2=2.*(th(iz+1)-th(iz))/(dz(iz+1)+dz(iz))                  
        dtmdz=0.5*(dtdz1+dtdz2)
        cdm=0.5*(exch(iz+1)+exch(iz))
        bu(iz)=-cdm*dtmdz*g/th0(iz)         
       enddo
!
                 
       bu(kte)=0.
         
       return
       end subroutine buoy

! ===6=8===============================================================72
! ===6=8===============================================================72

       subroutine shear(ix,iy,g,kms,kme,kts,kte,u,v,cdua,dz,sh,ustar,tstar,th,ufrac_int)
! compute shear term
       implicit none
       integer kms,kme,kts,kte,iz,ix,iy
       real dudz1,dudz2,dvdz1,dvdz2,cdm,dumdz,ustar
       real tstar,th,al,phim,g      
       real u(kms:kme),v(kms:kme),cdua(kms:kme),dz(kms:kme),sh(kms:kme)
       real u1,u2,v1,v2,ufrac_int

!       al=vk*g*tstar/(th*(ustar**2.))
!       if(al.ge.0.)phim=1.+4.7*dz(1)/2.*al
!       if(al.lt.0.)phim=(1.-15*dz(1)/2.*al)**(-0.25)       
!        
!       sh(1)=(ustar**3.)/vk/(dz(1)/2.)*(1.-ufrac_int)       
       sh(kts)=0.
       do iz=kts+1,kte-1
        u2=(dz(iz+1)*u(iz)+dz(iz)*u(iz+1))/(dz(iz)+dz(iz+1))
        u1=(dz(iz)*u(iz-1)+dz(iz-1)*u(iz))/(dz(iz-1)+dz(iz))
        v2=(dz(iz+1)*v(iz)+dz(iz)*v(iz+1))/(dz(iz)+dz(iz+1))
        v1=(dz(iz)*v(iz-1)+dz(iz-1)*v(iz))/(dz(iz-1)+dz(iz))        
        cdm=0.5*(cdua(iz)+cdua(iz+1)) 
        dumdz=((u2-u1)/dz(iz))**2.+((v2-v1)/dz(iz))**2.            
        sh(iz)=cdm*dumdz                          
       enddo

!!!!!!!
       sh(kte)=0.
       
       return
       end subroutine shear

! ===6=8===============================================================72
! ===6=8===============================================================72

       subroutine invert(kms,kme,kts,kte,a,c,x)
       
!ccccccccccccccccccccccccccccccc       
! Aim: Inversion and resolution of a tridiagonal matrix
!          A X = C
! Input:
!  a(*,1) lower diagonal (Ai,i-1)
!  a(*,2) principal diagonal (Ai,i)
!  a(*,3) upper diagonal (Ai,i+1)
!  c      
! Output
!  x     results
!ccccccccccccccccccccccccccccccc

       implicit none
       integer in
       integer kts,kte,kms,kme
       real a(kms:kme,3),c(kms:kme),x(kms:kme)                       
        
        do in=kte-1,kts,-1                 
         c(in)=c(in)-a(in,3)*c(in+1)/a(in+1,2)
         a(in,2)=a(in,2)-a(in,3)*a(in+1,1)/a(in+1,2)
        enddo
        
        do in=kts+1,kte        
         c(in)=c(in)-a(in,1)*c(in-1)/a(in-1,2)
        enddo
        
        do in=kts,kte
          
         x(in)=c(in)/a(in,2)
          
        enddo

        return
        end subroutine invert
  
! ===6=8===============================================================72
! ===6=8===============================================================72
              
       subroutine tke_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,vl,u,v,th,te,th0,ustar,tstar,exch,   &
                         dls,td,sh,bu,b_e,a_e,sf,ufrac_int)
! in this routine the shear, buoyancy and part of the dissipation terms
! of the TKE equation are computed       

       implicit none
       integer kms,kme,kts,kte,iz,ix,iy
       real g,ustar,tstar,ufrac_int
       real z(kms:kme),dz(kms:kme),u(kms:kme),v(kms:kme),th(kms:kme),th0(kms:kme),te(kms:kme)
       real exch(kms:kme),dls(kms:kme),td(kms:kme),sh(kms:kme),bu(kms:kme)
       real a_e(kms:kme),b_e(kms:kme)
       real vl(kms:kme),sf(kms:kme)
       real te1,dl1
    
       call shear(ix,iy,g,kms,kme,kts,kte,u,v,exch,dz,sh,ustar,tstar,th(kts),ufrac_int)
      
       call buoy(ix,iy,g,kms,kme,kts,kte,th,th0,exch,dz,bu,ustar,tstar,ufrac_int)
     
       do iz=kts,kte          
        te1=max(te(iz),temin)
        dl1=max(dls(iz),0.1)
        td(iz)=-ceps_b*sqrt(te1)/dl1
        sh(iz)=sh(iz)*sf(iz)
        bu(iz)=bu(iz)*sf(iz)
        a_e(iz)=a_e(iz)+td(iz)       
        b_e(iz)=b_e(iz)+sh(iz)+bu(iz)
       enddo 

         
       return
       end subroutine tke_bougeault    
!######################################################################
       subroutine pbl_height(kms,kme,kts,kte,dz,z,th,q,pblh)

! this routine computes the PBL height
! with an approach similar to MYNN
       implicit none
       integer kms,kme,kts,kte,iz
       real z(kms:kme),dz(kms:kme),th(kms:kme),q(kms:kme)
       real pblh
!Local
       real thv(kms:kme),zc(kms:kme)
       real thsfc

! compute the height of the center of the grid cells
      do iz=kts,kte
       zc(iz)=z(iz)+dz(iz)/2.
      enddo

! compute the virtual potential temperature

       do iz=kts,kte
        thv(iz)=th(iz)*(1.+0.61*q(iz))
       enddo
! now compute the PBL height

       pblh=0.
       thsfc=thv(kts)+0.5
       do iz=kts+1,kte
        if(pblh.eq.0.and.thv(iz).gt.thsfc)then
         pblh=zc(iz-1)+(thsfc-thv(iz-1))/(max(0.01,thv(iz)-thv(iz-1)))*(zc(iz)-zc(iz-1))
!         pblh=z(iz-1)+(thsfc-thv(iz-1))/(max(0.01,thv(iz)-thv(iz-1)))*(z(iz)-z(iz-1))
        endif
       enddo

       return
       end subroutine pbl_height


! ===6=8===============================================================72


! ===6=8===============================================================72
      SUBROUTINE BOULACINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,          &
     &                      TKE_PBL,EXCH_H,RESTART,ALLOWED_TO_READ,     &
     &                      IDS,IDE,JDS,JDE,KDS,KDE,                    &
     &                      IMS,IME,JMS,JME,KMS,KME,                    &
     &                      ITS,ITE,JTS,JTE,KTS,KTE                 )
!-----------------------------------------------------------------------
      IMPLICIT NONE
!-----------------------------------------------------------------------
      LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART
      INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE,                    &
     &                      IMS,IME,JMS,JME,KMS,KME,                    &
     &                      ITS,ITE,JTS,JTE,KTS,KTE

      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) ::    EXCH_H, &
     &                                                         RUBLTEN, &
     &                                                         RVBLTEN, &
     &                                                        RTHBLTEN, &
     &                                                        RQVBLTEN, &
     &                                                         TKE_PBL
      INTEGER :: I,J,K,ITF,JTF,KTF
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------

      JTF=MIN0(JTE,JDE-1)
      KTF=MIN0(KTE,KDE-1)
      ITF=MIN0(ITE,IDE-1)

      IF(.NOT.RESTART)THEN
        DO J=JTS,JTF
        DO K=KTS,KTF
        DO I=ITS,ITF
          TKE_PBL(I,K,J)=0.0001
          RUBLTEN(I,K,J)=0.
          RVBLTEN(I,K,J)=0.
          RTHBLTEN(I,K,J)=0.
          RQVBLTEN(I,K,J)=0.
          EXCH_H(I,K,J)=0.
        ENDDO
        ENDDO
        ENDDO
      ENDIF

      END SUBROUTINE BOULACINIT



END MODULE module_bl_boulac
Back to Top