-
Notifications
You must be signed in to change notification settings - Fork 38
/
ccsm_flux.F
1977 lines (1868 loc) · 78.9 KB
/
ccsm_flux.F
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
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#include "cppdefs.h"
MODULE ccsm_flux_mod
#if defined BULK_FLUXES || defined BULK_FLUXES2D
#if defined CCSM_FLUXES
!
!svn $Id$
!================================================== Hernan G. Arango ===
! Copyright (c) 2002-2020 The ROMS/TOMS Group !
! Licensed under a MIT/X style license !
! See License_ROMS.txt !
!=======================================================================
! !
! This routine computes the bulk parameterization of surface wind !
! stress and surface net heat fluxes. !
! !
! References: !
! !
! CCSM's flux_mod.F90 by B. Kauffman !
! !
!=======================================================================
!
implicit none
PRIVATE
! !PUBLIC MEMBER FUNCTIONS:
PUBLIC :: ccsm_flux
contains
!
!***********************************************************************
SUBROUTINE ccsm_flux (ng, tile)
!***********************************************************************
!
USE mod_param
USE mod_forces
USE mod_grid
USE mod_mixing
USE mod_ocean
# ifdef ICE_MODEL
USE mod_ice
# endif
# ifdef EVAP_SHALLOW
USE mod_coupling
# endif
# if defined BEST_NPZ && defined CLIM_ICE_1D
USE mod_clima
# endif
USE mod_stepping
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile
!
! Local variable declarations.
!
# include "tile.h"
!
# ifdef PROFILE
CALL wclock_on (ng, iNLM, 17, __LINE__, __FILE__)
# endif
CALL ccsm_flux_tile (ng, tile, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
& nrhs(ng), &
# ifdef ICE_MODEL
& liold(ng), &
& linew(ng), &
# endif
# if defined BEST_NPZ && defined ICE_BIO
& nstp(ng), &
# endif
# ifdef MASKING
& GRID(ng) % rmask, &
& GRID(ng) % umask, &
& GRID(ng) % vmask, &
# endif
# ifdef WET_DRY
& GRID(ng) % rmask_wet, &
& GRID(ng) % umask_wet, &
& GRID(ng) % vmask_wet, &
# endif
& OCEAN(ng) % t, &
# ifdef WIND_MINUS_CURRENT
& OCEAN(ng) % u, &
& OCEAN(ng) % v, &
# endif
& FORCES(ng) % Hair, &
& FORCES(ng) % Pair, &
& FORCES(ng) % Tair, &
& FORCES(ng) % Uwind, &
& FORCES(ng) % Vwind, &
# ifdef CLOUDS
& FORCES(ng) % cloud, &
# endif
# ifdef COARE_TAYLOR_YELLAND
& FORCES(ng) % Hwave, &
# endif
# if defined COARE_TAYLOR_YELLAND || defined COARE_OOST
& FORCES(ng) % Pwave_top, &
# endif
# if !defined DEEPWATER_WAVES && \
(defined COARE_TAYLOR_YELLAND || defined COARE_OOST)
& FORCES(ng) % Lwave, &
# endif
# ifdef RUNOFF
& FORCES(ng) % runoff, &
# endif
# ifdef EVAP_SHALLOW
& GRID(ng) % h, &
& COUPLING(ng) % Zt_avg1, &
# endif
& FORCES(ng) % rain, &
& FORCES(ng) % lhflx, &
& FORCES(ng) % lrflx, &
& FORCES(ng) % shflx, &
& FORCES(ng) % srflx, &
& FORCES(ng) % snhflx, &
# ifdef CICE_MODEL
& FORCES(ng) % LW_down, &
& FORCES(ng) % SW_down, &
# endif
# if defined ALBEDO && defined SHORTWAVE
& FORCES(ng) % albedo, &
# ifdef ICE_MODEL
& FORCES(ng) % albedo_ice, &
# ifdef CCSM_ICE_SHORTWAVE
& FORCES(ng) % frswpen, &
# endif
# endif
# endif
# if defined SHORTWAVE && defined ALBEDO_DIRDIFF
& FORCES(ng) % cawdir, &
# endif
# ifdef ICE_MODEL
& ICE(ng) % ai, &
& ICE(ng) % hi, &
& ICE(ng) % rfaci, & !!!Inserted PWA 25/09/2019
# ifdef ICE_THERMO
& FORCES(ng) % qai_n, &
& FORCES(ng) % qi_o_n, &
& FORCES(ng) % SW_thru_ice, &
& FORCES(ng) % qswi_n, &
& ICE(ng) % hsn, &
& ICE(ng) % tis, &
& ICE(ng) % coef_ice_heat, &
& ICE(ng) % rhs_ice_heat, &
& FORCES(ng) % snow_n, &
# ifdef SNOWFALL
& FORCES(ng) % snow, &
# endif
# endif
& FORCES(ng) % sustr_aw, &
& FORCES(ng) % svstr_aw, &
& FORCES(ng) % qao_n, &
& FORCES(ng) % tau_aix_n, &
& FORCES(ng) % tau_aiy_n, &
# endif
# if defined BEST_NPZ
# if defined CLIM_ICE_1D
& OCEAN(ng) % it, &
& CLIMA(ng) % tclm, &
# elif defined ICE_BIO
& OCEAN(ng) % it, &
& ICE(ng) % IcePhL, &
& ICE(ng) % IceLog, &
# endif
# endif
# ifdef EMINUSP
& FORCES(ng) % EminusP, &
& FORCES(ng) % evap, &
# endif
# ifdef ICE_DIAGS
& FORCES(ng) % LW_down, &
& FORCES(ng) % SW_down, &
& FORCES(ng) % lat_ice, &
& FORCES(ng) % sens_ice, &
& FORCES(ng) % LW_up_ice, &
& FORCES(ng) % SW_up_ice, &
& FORCES(ng) % saltflux_ocean, &
# endif
& FORCES(ng) % sustr, &
& FORCES(ng) % svstr)
# ifdef PROFILE
CALL wclock_off (ng, iNLM, 17, __LINE__, __FILE__)
# endif
RETURN
END SUBROUTINE ccsm_flux
!
!***********************************************************************
SUBROUTINE ccsm_flux_tile (ng, tile, &
& LBi, UBi, LBj, UBj, &
& IminS, ImaxS, JminS, JmaxS, &
& nrhs, &
# ifdef ICE_MODEL
& liold, linew, &
# endif
# if defined BEST_NPZ && defined ICE_BIO
& nstp, &
# endif
# ifdef MASKING
& rmask, umask, vmask, &
# endif
# ifdef WET_DRY
& rmask_wet, umask_wet, vmask_wet, &
# endif
& t, &
# ifdef WIND_MINUS_CURRENT
& u, v, &
# endif
& Hair, Pair, Tair, Uwind, Vwind, &
# ifdef CLOUDS
& cloud, &
# endif
# ifdef COARE_TAYLOR_YELLAND
& Hwave, &
# endif
# if defined COARE_TAYLOR_YELLAND || defined COARE_OOST
& Pwave_top, &
# endif
# if !defined DEEPWATER_WAVES && \
(defined COARE_TAYLOR_YELLAND || defined COARE_OOST)
& Lwave, &
# endif
# ifdef RUNOFF
& runoff, &
# endif
# ifdef EVAP_SHALLOW
& h, Zt_avg1, &
# endif
& rain, lhflx, lrflx, shflx, &
& srflx, snhflx, &
# ifdef CICE_MODEL
& LW_down, SW_down, &
# endif
# if defined ALBEDO && defined SHORTWAVE
& albedo, &
# ifdef ICE_MODEL
& albedo_ice, &
# ifdef CCSM_ICE_SHORTWAVE
& frswpen, &
# endif
# endif
# endif
# if defined SHORTWAVE && defined ALBEDO_DIRDIFF
& cawdir, &
# endif
# ifdef ICE_MODEL
& ai, hi, rfaci, & !!!Inserted PWA 25/09/2019, rfaci
# ifdef ICE_THERMO
& qai_n, qi_o_n, SW_thru_ice, qswi_n, &
& hsn, tis, coef_ice_heat, &
& rhs_ice_heat, snow_n, &
# ifdef SNOWFALL
& snow, &
# endif
# endif
& sustr_aw, svstr_aw, qao_n, &
& tau_aix_n, tau_aiy_n, &
# endif
# if defined BEST_NPZ
# if defined CLIM_ICE_1D
& it, tclm, &
# elif defined ICE_BIO
& it, IcePhL, IceLog, &
# endif
# endif
# ifdef EMINUSP
& EminusP, evap, &
# endif
# ifdef ICE_DIAGS
& LW_down, SW_down, &
& lat_ice, sens_ice, &
& LW_up_ice, SW_up_ice, &
& saltflux_ocean, &
# endif
& sustr, svstr)
!***********************************************************************
! call srfflx_ao( nloc , z , u , v ,ptem , &
! & shum ,dens ,uocn ,vocn , &
! & tocn ,mask , shflx , lhflx ,lwup , &
! & evap ,taux , tauy ,tref ,qref , &
! & duu10n )
!***********************************************************************
!
USE mod_param
USE mod_scalars
# if defined BEST_NPZ
USE mod_biology
# endif
# if defined BEST_NPZ && defined CLIM_ICE_1D
USE mod_clima
# endif
# ifdef ICE_BOX
USE dateclock_mod, ONLY : caldate
# endif
USE exchange_2d_mod
# ifdef DISTRIBUTE
USE mp_exchange_mod, ONLY : mp_exchange2d
# endif
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile
integer, intent(in) :: LBi, UBi, LBj, UBj
integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
integer, intent(in) :: nrhs
# ifdef ICE_MODEL
integer, intent(in) :: liold
integer, intent(in) :: linew
# endif
# if defined BEST_NPZ && defined ICE_BIO
integer, intent(in) :: nstp
# endif
!
# ifdef ASSUMED_SHAPE
# ifdef MASKING
real(r8), intent(in) :: rmask(LBi:,LBj:)
real(r8), intent(in) :: umask(LBi:,LBj:)
real(r8), intent(in) :: vmask(LBi:,LBj:)
# endif
# ifdef WET_DRY
real(r8), intent(in) :: rmask_wet(LBi:,LBj:)
real(r8), intent(in) :: umask_wet(LBi:,LBj:)
real(r8), intent(in) :: vmask_wet(LBi:,LBj:)
# endif
real(r8), intent(in) :: t(LBi:,LBj:,:,:,:)
# ifdef WIND_MINUS_CURRENT
real(r8), intent(in) :: u(LBi:,LBj:,:,:)
real(r8), intent(in) :: v(LBi:,LBj:,:,:)
# endif
real(r8), intent(in) :: Hair(LBi:,LBj:)
real(r8), intent(in) :: Pair(LBi:,LBj:)
real(r8), intent(in) :: Tair(LBi:,LBj:)
# ifdef WIND_MINUS_CURRENT
real(r8), intent(inout) :: Uwind(LBi:,LBj:)
real(r8), intent(inout) :: Vwind(LBi:,LBj:)
# else
real(r8), intent(in) :: Uwind(LBi:,LBj:)
real(r8), intent(in) :: Vwind(LBi:,LBj:)
# endif
# ifdef CLOUDS
real(r8), intent(in) :: cloud(LBi:,LBj:)
# endif
# ifdef COARE_TAYLOR_YELLAND
real(r8), intent(in) :: Hwave(LBi:,LBj:)
# endif
# if defined COARE_TAYLOR_YELLAND || defined COARE_OOST
real(r8), intent(in) :: Pwave_top(LBi:,LBj:)
# endif
# if !defined DEEPWATER_WAVES && \
(defined COARE_TAYLOR_YELLAND || defined COARE_OOST)
real(r8), intent(in) :: Lwave(LBi:,LBj:)
# endif
real(r8), intent(in) :: rain(LBi:,LBj:)
# ifdef RUNOFF
real(r8), intent(in) :: runoff(LBi:,LBj:)
# endif
# ifdef EVAP_SHALLOW
real(r8), intent(in) :: h(LBi:,LBj:)
real(r8), intent(in) :: Zt_avg1(LBi:,LBj:)
# endif
real(r8), intent(inout) :: lhflx(LBi:,LBj:)
real(r8), intent(inout) :: lrflx(LBi:,LBj:)
real(r8), intent(inout) :: shflx(LBi:,LBj:)
real(r8), intent(inout) :: srflx(LBi:,LBj:)
real(r8), intent(inout) :: snhflx(LBi:,LBj:)
# ifdef CICE_MODEL
real(r8), intent(out) :: LW_down(LBi:,LBj:)
real(r8), intent(out) :: SW_down(LBi:,LBj:)
# endif
# if defined ALBEDO && defined SHORTWAVE
real(r8), intent(inout) :: albedo(LBi:,LBj:)
# ifdef ICE_MODEL
real(r8), intent(inout) :: albedo_ice(LBi:,LBj:)
# ifdef CCSM_ICE_SHORTWAVE
real(r8), intent(in) :: frswpen(LBi:,LBj:)
# endif
# endif
# endif
# if defined SHORTWAVE && defined ALBEDO_DIRDIFF
real(r8), intent(inout) :: cawdir(LBi:,LBj:)
# endif
# ifdef ICE_MODEL
real(r8), intent(out) :: sustr_aw(LBi:,LBj:)
real(r8), intent(out) :: svstr_aw(LBi:,LBj:)
real(r8), intent(out) :: qao_n(LBi:,LBj:)
real(r8), intent(in) :: ai(LBi:,LBj:,:)
real(r8), intent(in) :: hi(LBi:,LBj:,:)
real(r8), intent(out) :: rfaci(LBi:,LBj:)
# ifdef ICE_THERMO
real(r8), intent(out) :: qai_n(LBi:,LBj:)
real(r8), intent(out) :: qi_o_n(LBi:,LBj:)
real(r8), intent(out) :: SW_thru_ice(LBi:,LBj:)
real(r8), intent(out) :: qswi_n(LBi:,LBj:)
real(r8), intent(in) :: hsn(LBi:,LBj:,:)
real(r8), intent(in) :: tis(LBi:,LBj:)
real(r8), intent(out) :: coef_ice_heat(LBi:,LBj:)
real(r8), intent(out) :: rhs_ice_heat(LBi:,LBj:)
real(r8), intent(out) :: snow_n(LBi:,LBj:)
# ifdef SNOWFALL
real(r8), intent(in) :: snow(LBi:,LBj:)
# endif
# endif
real(r8), intent(out) :: tau_aix_n(LBi:,LBj:)
real(r8), intent(out) :: tau_aiy_n(LBi:,LBj:)
# endif
# if defined BEST_NPZ
# if defined CLIM_ICE_1D
real(r8), intent(in) ::tclm(LBi:,LBj:,:,:,:)
real(r8), intent(in) :: it(LBi:,LBj:,:,:)
# elif defined ICE_BIO
real(r8), intent(in) :: it(LBi:,LBj:,:,:)
real(r8), intent(in) :: IcePhL(LBi:,LBj:,:)
integer, intent(inout) :: IceLog(LBi:,LBj:,:)
# endif
# endif
# ifdef EMINUSP
real(r8), intent(out) :: EminusP(LBi:,LBj:)
real(r8), intent(out) :: evap(LBi:,LBj:)
# endif
# ifdef ICE_DIAGS
real(r8), intent(out) :: LW_down(LBi:,LBj:)
real(r8), intent(out) :: SW_down(LBi:,LBj:)
real(r8), intent(out) :: lat_ice(LBi:,LBj:)
real(r8), intent(out) :: sens_ice(LBi:,LBj:)
real(r8), intent(out) :: LW_up_ice(LBi:,LBj:)
real(r8), intent(out) :: SW_up_ice(LBi:,LBj:)
real(r8), intent(out) :: saltflux_ocean(LBi:,LBj:)
# endif
real(r8), intent(out) :: sustr(LBi:,LBj:)
real(r8), intent(out) :: svstr(LBi:,LBj:)
# else
# ifdef MASKING
real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
# endif
# ifdef WET_DRY
real(r8), intent(in) :: rmask_wet(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: umask_wet(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: vmask_wet(LBi:UBi,LBj:UBj)
# endif
real(r8), intent(in) :: t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
# ifdef WIND_MINUS_CURRENT
real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),3)
real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),3)
# endif
real(r8), intent(in) :: Hair(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: Pair(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: Tair(LBi:UBi,LBj:UBj)
# ifdef WIND_MINUS_CURRENT
real(r8), intent(inout) :: Uwind(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: Vwind(LBi:UBi,LBj:UBj)
# else
real(r8), intent(in) :: Uwind(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: Vwind(LBi:UBi,LBj:UBj)
# endif
# ifdef CLOUDS
real(r8), intent(in) :: cloud(LBi:UBi,LBj:UBj)
# endif
# ifdef COARE_TAYLOR_YELLAND
real(r8), intent(in) :: Hwave(LBi:UBi,LBj:UBj)
# endif
# if defined COARE_TAYLOR_YELLAND || defined COARE_OOST
real(r8), intent(in) :: Pwave_top(LBi:UBi,LBj:UBj)
# endif
# if !defined DEEPWATER_WAVES && \
(defined COARE_TAYLOR_YELLAND || defined COARE_OOST)
real(r8), intent(in) :: Lwave(LBi:UBi,LBj:UBj)
# endif
real(r8), intent(in) :: rain(LBi:UBi,LBj:UBj)
# ifdef RUNOFF
real(r8), intent(in) :: runoff(LBi:UBi,LBj:UBj)
# endif
# ifdef EVAP_SHALLOW
real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: Zt_avg1(LBi:UBi,LBj:UBj)
# endif
real(r8), intent(inout) :: lhflx(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: lrflx(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: shflx(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: srflx(LBi:UBi,LBj:UBj)
real(r8), intent(inout) :: snhflx(LBi:UBi,LBj:UBj)
# ifdef CICE_MODEL
real(r8), intent(out) :: LW_down(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: SW_down(LBi:UBi,LBj:UBj)
# endif
# if defined ALBEDO && defined SHORTWAVE
real(r8), intent(inout) :: albedo(LBi:UBi,LBj:UBj)
# ifdef ICE_MODEL
real(r8), intent(inout) :: albedo_ice(LBi:UBi,LBj:UBj)
# ifdef CCSM_ICE_SHORTWAVE
real(r8), intent(in) :: frswpen(LBi:UBi,LBj:UBj)
# endif
# endif
# endif
# if defined SHORTWAVE && defined ALBEDO_DIRDIFF
real(r8), intent(inout) :: cawdir(LBi:UBi,LBj:UBj)
# endif
# ifdef ICE_MODEL
real(r8), intent(out) :: sustr_aw(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: svstr_aw(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: qao_n(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: ai(LBi:UBi,LBj:UBj,2)
real(r8), intent(in) :: hi(LBi:UBi,LBj:UBj,2)
real(r8), intent(out) :: rfaci(LBi:UBi,LBj:UBj)
# ifdef ICE_THERMO
real(r8), intent(out) :: qai_n(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: qi_o_n(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: SW_thru_ice(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: qswi_n(LBi:UBi,LBj:UBj)
real(r8), intent(in) :: hsn(LBi:UBi,LBj:UBj,2)
real(r8), intent(in) :: tis(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: coef_ice_heat(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: rhs_ice_heat(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: snow_n(LBi:UBi,LBj:UBj)
# ifdef SNOWFALL
real(r8), intent(in) :: snow(LBi:UBi,LBj:UBj)
# endif
# endif
real(r8), intent(out) :: tau_aix_n(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: tau_aiy_n(LBi:UBi,LBj:UBj)
# endif
# if defined BEST_NPZ
# if defined CLIM_ICE_1D
real(r8), intent(in) :: it(LBi:UBi,LBj:UBj,3,1)
real(r8), intent(in) ::tclm(LBi:UBi,LBj:UBj,UBk,3,UBt+1)
# elif defined ICE_BIO
real(r8), intent(in) :: it(LBi:UBi,LBj:UBj,3,1)
real(r8), intent(in) :: IcePhL(LBi:,LBj:,:)
integer, intent(inout) :: IceLog(LBi:,LBj:,:)
# endif
# endif
# ifdef EMINUSP
real(r8), intent(out) :: EminusP(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: evap(LBi:UBi,LBj:UBj)
# endif
# ifdef ICE_DIAGS
real(r8), intent(out) :: LW_down(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: SW_down(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: lat_ice(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: sens_ice(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: LW_up_ice(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: SW_up_ice(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: saltflux_ocean(LBi:UBi,LBj:UBj)
# endif
real(r8), intent(out) :: sustr(LBi:UBi,LBj:UBj)
real(r8), intent(out) :: svstr(LBi:UBi,LBj:UBj)
# endif
!
! Local variable declarations.
!
integer :: i, j, listp
integer :: month
real(r8), parameter :: eps = 1.0E-20_r8
real(r8), parameter :: r3 = 1.0_r8/3.0_r8
real(r8), parameter :: Ch = 1.75e-3
real(r8), parameter :: Ce = 1.75e-3
! ratio of molecular weight of water to dry air
real(r8), parameter :: epsilon = 0.622_r8
real(r8) :: Hscale, Hscale2
real(r8) :: RH, cff1, cff2
real(r8) :: cff
# ifdef ICE_MODEL
real(r8) :: transw, transi, rtransi
# endif
# ifdef ICE_BULK_FLUXES
real(r8) :: qswi, qlwi, qlh_i, qsh_i
real(r8) :: le_i, dq_i,fqlat1, sfc_temp, slp, Qsati
real(r8) :: vap_p_i
real(r8), parameter :: Io_frac = 0.17_r8 ! from MU
# endif
# if defined SHORTWAVE && defined ALBEDO_DIRDIFF
real(r8) :: cawdif
real(r8), parameter :: albw_d=0.065_r8
# endif
# ifdef EVAP_SHALLOW
real(r8), parameter :: coef_shallow = 5.28_r8
! real(r8), parameter :: coef_shallow = 4.69_r8
real(r8) :: cffe
# endif
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Uair
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vair
real(r8) :: fsnow
# ifdef ICE_BOX
real(r8), parameter :: sensible_MU(12) = &
& (/ 19.1, 12.3, 11.6, 4.7, -7.3, -6.3, &
& -4.8, -6.5, -2.7, 1.6, 9.0, 12.8 /)
real(r8), parameter :: latent_MU(12) = &
& (/ 0.0, -0.32, -0.48, -1.45, -7.4, -11.3, &
& -10.3, -10.7, -6.3, -3.1, -.16, -.16 /)
# endif
! ------------------------------------------------------------------
! !DESCRIPTION:
! wrapper to atm/ocn flux calculation
!
! !REMARKS:
! All data must be on the ocean domain (note: a domain includes a
! particular decomposition).
!
! !REVISION HISTORY:
! 2002-Jun-10 - B. Kauffman - first version
!
! ------------------------------------------------------------------
!
! !IROUTINE: srfflx_ao -- internal atm/ocn flux calculation
!
! !DESCRIPTION:
!
! Internal atm/ocn flux calculation
!
! !REVISION HISTORY:
! 2002-Jun-10 - B. Kauffman - brought in from cpl5.
! 2003-Apr-02 - B. Kauffman - taux & tauy now utilize ocn velocity
! 2003-Apr-02 - B. Kauffman - tref,qref,duu10n mods as per Bill Large
!
! !INTERFACE: ------------------------------------------------------------------
!SUBROUTINE srfflx_ao_tile(blk_ZX ,Uwind ,Vwind ,Tair , &
! & Hair ,rhoAir ,u ,v , &
! & t ,rmask ,shflx ,lhflx ,lwup , &
! & evap ,taux ,tauy ,tref ,qref , &
! & duu10n )
! !INPUT/OUTPUT PARAMETERS:
!--- input arguments --------------------------------
! integer(IN),intent(in) :: rmask(imax) ! ocn domain mask 0
! real(R8) ,intent(in) :: blk_ZX (imax) ! atm level height (m)
! real(R8) ,intent(in) :: Uwind (imax) ! atm u wind (m/s)
! real(R8) ,intent(in) :: Vwind (imax) ! atm v wind (m/s)
! real(R8) ,intent(in) :: Tair(imax) ! atm potential T (K)
! real(R8) ,intent(in) :: Hair (imax) ! atm specific humidity (kg/kg)
! real(R8) ,intent(in) :: rhoAir (imax) ! atm air density (kg/m^3)
! real(R8) ,intent(in) :: u (imax) ! ocn u-velocity (m/s)
! real(R8) ,intent(in) :: v (imax) ! ocn v-velocity (m/s)
! real(R8) ,intent(in) :: t (imax) ! ocn temperature (K)
!--- output arguments -------------------------------
! real(R8),intent(out) :: shflx (imax) ! heat flux: sensible (W/m^2)
! real(R8),intent(out) :: lhflx (imax) ! heat flux: latent (W/m^2)
! real(R8),intent(out) :: lwup (imax) ! heat flux: lw upward (W/m^2)
! real(R8),intent(out) :: evap (imax) ! water flux: evap ((kg/s)/m^2)
! real(R8),intent(out) :: taux (imax) ! surface stress, zonal (N)
! real(R8),intent(out) :: tauy (imax) ! surface stress, maridional (N)
! real(R8),intent(out) :: tref (imax) ! diag: 2m ref height T (K)
! real(R8),intent(out) :: qref (imax) ! diag: 2m ref humidity (kg/kg)
! real(R8),intent(out) :: duu10n(imax) ! diag: 10m wind speed squared (m/s)^2
!EOP
!--- local constants --------------------------------
real(R8),parameter :: umin = 0.5 ! minimum wind speed (m/s)
real(R8),parameter :: zref = 10.0 ! reference height (m)
real(R8),parameter :: ztref = 2.0 ! reference height for air T (m)
!--- local variables --------------------------------
real(r8) :: lwup ! upward longwave radiation
real(r8) :: vmag ! surface wind magnitude (m/s)
real(r8) :: thvbot ! virtual temperature (K)
real(r8) :: ssq ! sea surface humidity (kg/kg)
real(r8) :: Tseak ! SST in Kelvins (K)
real(r8) :: delt ! potential T difference (K)
real(r8) :: delq ! humidity difference (kg/kg)
real(r8) :: stable ! stability factor
real(r8) :: rdn ! sqrt of neutral exchange coeff (momentum)
real(r8) :: rhn ! sqrt of neutral exchange coeff (heat)
real(r8) :: ren ! sqrt of neutral exchange coeff (water)
real(r8) :: rd ! sqrt of exchange coefficient (momentum)
real(r8) :: re ! sqrt of exchange coefficient (water)
real(r8) :: ustar ! ustar
real(r8) :: qstar ! qstar
real(r8) :: tstar ! tstar
real(r8) :: hol ! H (at blk_ZX) over L
real(r8) :: xsq ! ?
real(r8) :: xqq ! ?
real(r8) :: psimh ! stability function at blk_ZX (momentum)
real(r8) :: psixh ! stability function at blk_ZX (heat and water)
real(r8) :: psix2 ! stability function at ztref reference height
real(r8) :: alz ! ln(blk_ZX/zref)
real(r8) :: al2 ! ln(zref/ztref)
real(r8) :: u10n ! 10m neutral wind
real(r8) :: tau ! stress at blk_ZX
real(r8) :: cpair ! specific heat of moist air
real(r8) :: fac ! vertical interpolation factor
real(r8) :: Hlv ! Latent heat of vaporization at sea surface
real(r8), dimension(IminS:ImaxS) :: rhoAir
real(r8), dimension(IminS:ImaxS) :: TairK
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: LHeat
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: SHeat
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Taux
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Tauy
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: u10n_tmp
# ifdef ICE_THERMO
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Coef_Ice
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: RHS_Ice
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Qai
# endif
# ifdef ICE_MODEL
real(r8), parameter :: fv=0.43_r8 !visual fraction of total shortwave (Aumont et al., 2015)
# ifdef CCSM_ICE_SHORTWAVE
real(r8), parameter :: kappav = 1.4_r8 !vis extnctn coef in ice (1/m); CESM2/CICE5.1 default, DuVivier et al., 2018 (cf. 0.73 in Briegleb et al., 2004)
# endif
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Taux_Ice
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Tauy_Ice
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ice_thick
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: snow_thick
real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: i0frac
# endif
!--- local functions --------------------------------
real(r8) :: qsat ! function: the saturation humididty of air (kg/m^3)
real(r8) :: cdn ! function: neutral drag coeff at 10m
real(r8) :: psimhu ! function: unstable part of psimh
real(r8) :: psixhu ! function: unstable part of psimx
real(r8) :: tsk ! water temperature (K)
real(r8) :: Tk ! dummy arg ~ air temperature (K)
real(r8) :: Umps ! dummy arg ~ wind velocity (m/s)
real(r8) :: xd ! dummy arg ~ ?
! Stefan-Boltzmann constant ~ W/m^2/K^4
real(r8),parameter :: shr_const_stebol = 5.67e-8_r8
! Boltzmann's constant ~ J/K/molecule
real(r8),parameter :: shr_const_boltz = 1.38065e-23_r8
! Avogadro's number ~ molecules/kmole
real(r8),parameter :: shr_const_avogad = 6.02214e26_r8
! Universal gas constant ~ J/K/kmole
real(r8),parameter :: shr_const_rgas = &
& shr_const_avogad*shr_const_boltz
! molecular weight dry air ~ kg/kmole
real(r8),parameter :: shr_const_mwdair = 28.966_r8
! molecular weight water vapor
real(r8),parameter :: shr_const_mwwv = 18.016_r8
! Water vapor gas constant ~ J/K/kg
real(r8),parameter :: shr_const_rwv = &
& shr_const_rgas/shr_const_mwwv
! Dry air gas constant ~ J/K/kg
real(r8),parameter :: shr_const_rdair = &
& shr_const_rgas/shr_const_mwdair
real(r8),parameter :: shr_const_zvir = &
& (shr_const_rwv/shr_const_rdair)-1.0_r8
! Von Karman constant
real(r8),parameter :: shr_const_karman = 0.4_r8
! specific heat of dry air ~ J/kg/K
real(r8),parameter :: shr_const_cpdair = 1.00464e3_r8
! specific heat of water vap ~ J/kg/K
real(r8),parameter :: shr_const_cpwv = 1.810e3_R8
! CPWV/CPDAIR - 1.0
real(r8),parameter :: shr_const_cpvir = &
& (shr_const_cpwv/shr_const_cpdair)-1.0_r8
! latent heat of evaporation ~ J/kg
real(r8),parameter :: shr_const_latvap = 2.501e6_r8
qsat(Tk) = 640380.0 / exp(5107.4/Tk)
cdn(Umps) = 0.0027 / Umps + 0.000142 + 0.0000764 * Umps
psimhu(xd) = log((1.0+xd*(2.0+xd))*(1.0+xd*xd)/8.0) - &
& 2.0*atan(xd) + 1.571
psixhu(xd) = 2.0 * log((1.0 + xd*xd)/2.0)
#include "set_bounds.h"
!-------------------------------------------------------------------------------
! PURPOSE:
! computes atm/ocn surface fluxes
!
! NOTES:
! o all fluxes are positive downward
! o net heat flux = net sw + lw up + lw down + shflx + lhflx
! o here, tstar = <WT>/U*, and qstar = <WQ>/U*.
! o wind speeds should all be above a minimum speed (eg. 1.0 m/s)
!
! ASSUMPTIONS:
! o Neutral 10m drag coeff: cdn = .0027/U10 + .000142 + .0000764 U10
! o Neutral 10m stanton number: ctn = .0327 sqrt(cdn), unstable
! ctn = .0180 sqrt(cdn), stable
! o Neutral 10m dalton number: cen = .0346 sqrt(cdn)
! o The saturation humidity of air at T(K): qsat(T) (kg/m^3)
!-------------------------------------------------------------------------------
al2 = log(zref/ztref)
# if defined ICE_MODEL
IF (PerfectRST(ng) .and. iic(ng).eq.ntstart(ng)) THEN
listp = liold
ELSE
listp = linew
END IF
# endif
# ifdef ICE_BOX
CALL caldate(tdays(ng), mmi=month)
# endif
# ifdef WIND_MINUS_CURRENT
!
! Modify near-surface (2m or 10m) effective winds by subtracting the
! ocean surface current (J Wilkin). See:
!
! Bye, J.A.T. and J.-O. Wolff, 1999: Atmosphere-ocean momentum exchange
! in general circulation models. J. Phys. Oceanogr. 29, 671-691.
!
DO j=Jstr-1,JendR
DO i=Istr-1,IendR
Uair(i,j)=Uwind(i,j)- &
& 0.5_r8*(u(i,j,N(ng),nrhs)+u(i+1,j,N(ng),nrhs))
Vair(i,j)=Vwind(i,j)- &
& 0.5_r8*(v(i,j,N(ng),nrhs)+v(i,j+1,N(ng),nrhs))
END DO
END DO
# else
!
! Load wind components to local arrays.
!
DO j=Jstr-1,JendR
DO i=Istr-1,IendR
Uair(i,j)=Uwind(i,j)
Vair(i,j)=Vwind(i,j)
END DO
END DO
# endif
Hscale = rho0*Cp
Hscale2 = 1.0_r8/(rho0*Cp)
! Note that this loop needs to be cleaned of all global arrays for
! OpenMP.
DO j=JstrR,JendR
DO i=IstrR,IendR
!--- compute some needed quantities ---
vmag = max(umin, sqrt( (Uair(i,j)**2 + Vair(i,j)**2)) )
TairK(i) = Tair(i,j) + 273.16_r8
TseaK = t(i,j,N(ng),nrhs,itemp) + 273.16_r8
thvbot = TairK(i) * (1.0 + shr_const_zvir * Hair(i,j)) ! virtual temp (K)
rhoAir(i) = Pair(i,j)*100.0_r8/(blk_Rgas*TairK(i)* &
& (1.0_r8+0.61_r8*Hair(i,j))) ! Moist air density (kg/m3).
ssq = 0.98_r8 * qsat(TseaK) / rhoAir(i) ! sea surf hum (kg/kg)
delt = Tair(i,j) - t(i,j,N(ng),nrhs,itemp) ! pot temp diff (K)
delq = Hair(i,j) - ssq ! spec hum dif (kg/kg)
alz = log(blk_ZW(ng)/zref)
cpair = shr_const_cpdair*(1.0 + shr_const_cpvir*ssq)
!------------------------------------------------------------
! first estimate of Z/L and ustar, tstar and qstar
!------------------------------------------------------------
!--- neutral coefficients, z/L = 0.0 ---
stable = 0.5_r8 + sign(0.5_r8 , delt)
rdn = sqrt(cdn(vmag))
rhn = (1.0_r8-stable) * 0.0327_r8 + stable * 0.018_r8
ren = 0.0346_r8
!--- ustar, tstar, qstar ---
ustar = rdn * vmag
tstar = rhn * delt
qstar = ren * delq
!--- compute stability & evaluate all stability functions ---
hol = shr_const_karman*g*blk_ZT(ng)* &
& (tstar/thvbot+qstar/(1.0_r8/shr_const_zvir+Hair(i,j)))/ &
& ustar**2
hol = sign( min(abs(hol),10.0_r8), hol )
stable = 0.5_r8 + sign(0.5_r8 , hol)
xsq = max(sqrt(abs(1.0_r8 - 16.0_r8*hol)) , 1.0_r8)
xqq = sqrt(xsq)
psimh = -5.0_r8*hol*stable + (1.0_r8-stable)*psimhu(xqq)
psixh = -5.0_r8*hol*stable + (1.0_r8-stable)*psixhu(xqq)
!--- shift wind speed using old coefficient ---
rd = rdn / (1.0_r8 + rdn/shr_const_karman*(alz-psimh))
u10n = vmag * rd / rdn
!--- update transfer coeffs at 10m and neutral stability ---
rdn = sqrt(cdn(u10n))
ren = 0.0346_r8
rhn = (1.0_r8-stable)*0.0327_r8 + stable * 0.018_r8
!--- shift all coeffs to measurement height and stability ---
rd = rdn / (1.0_r8 + rdn/shr_const_karman*(alz-psimh))
rh = rhn / (1.0_r8 + rhn/shr_const_karman*(alz-psixh))
re = ren / (1.0_r8 + ren/shr_const_karman*(alz-psixh))
!--- update ustar, tstar, qstar using updated, shifted coeffs --
ustar = rd * vmag
tstar = rh * delt
qstar = re * delq
!------------------------------------------------------------
! iterate to converge on Z/L, ustar, tstar and qstar
!------------------------------------------------------------
!--- compute stability & evaluate all stability functions ---
hol = shr_const_karman*g*blk_ZQ(ng)* &
& (tstar/thvbot+qstar/(1.0/shr_const_zvir+Hair(i,j)))/ &
& ustar**2
hol = sign( min(abs(hol),10.0_r8), hol )
stable = 0.5_r8 + sign(0.5_r8 , hol)
xsq = max(sqrt(abs(1.0_r8 - 16.0_r8*hol)) , 1.0_r8)
xqq = sqrt(xsq)
psimh = -5.0_r8*hol*stable + (1.0_r8-stable)*psimhu(xqq)
psixh = -5.0_r8*hol*stable + (1.0_r8-stable)*psixhu(xqq)
!--- shift wind speed using old coeffs ---
rd = rdn / (1.0_r8 + rdn/shr_const_karman*(alz-psimh))
u10n = vmag * rd/rdn
!--- update transfer coeffs at 10m and neutral stability ---
rdn = sqrt(cdn(u10n))
ren = 0.0346_r8
rhn = (1.0_r8 - stable)*0.0327_r8 + stable * 0.018_r8
!--- shift all coeffs to measurement height and stability ---
rd = rdn / (1.0_r8 + rdn/shr_const_karman*(alz-psimh))
rh = rhn / (1.0_r8 + rhn/shr_const_karman*(alz-psixh))
re = ren / (1.0_r8 + ren/shr_const_karman*(alz-psixh))
!--- update ustar, tstar, qstar using updated, shifted coeffs ---
ustar = rd * vmag
tstar = rh * delt
qstar = re * delq
!--- RD: save u10n to u10n_tmp of dim IminS:ImaxS,JminS:JmaxS
! otherwise there is an array overflow !!!
u10n_tmp(i,j) = u10n
!------------------------------------------------------------
! compute the fluxes
!------------------------------------------------------------
tau = rhoAir(i) * ustar * ustar
# if defined ICE_MODEL && defined ICE_THERMO
qswi_n(i,j) = srflx(i,j) * Hscale
# endif
!--- momentum flux ---
Taux(i,j) = tau * Uair(i,j) / vmag
Tauy(i,j) = tau * Vair(i,j) / vmag
# ifdef ICE_MODEL
Taux_Ice(i,j) = rhoAir(i)*cdai(ng)*Uwind(i,j)*vmag
Tauy_Ice(i,j) = rhoAir(i)*cdai(ng)*Vwind(i,j)*vmag
# ifdef ICE_THERMO
! Correct ocean net short-wave radiation for ice cover and convert to
! kinematic units
Coef_Ice(i,j) = 0._r8
RHS_Ice(i,j) = 0._r8
! Sensible heat flux over ice
IF (ai(i,j,listp).gt.min_a(ng)) THEN
sfc_temp = tis(i,j) + 273.16_r8
ELSE
sfc_temp = t(i,j,N(ng),nrhs,itemp) + 273.16_r8
ENDIF
# ifdef ICE_BOX
qsh_i = sensible_MU(month)
RHS_Ice(i,j) = RHS_Ice(i,j) + qsh_i
# else
qsh_i = rhoAir(i)*blk_Cpa*Ch*vmag* &
& (TairK(i)+0.0098*blk_ZT(ng)-sfc_temp)
Coef_Ice(i,j) = Coef_Ice(i,j) + &
& rhoAir(i)*blk_Cpa*Ch*vmag
RHS_Ice(i,j) = RHS_Ice(i,j) + &
& rhoAir(i)*blk_Cpa*Ch*vmag*(TairK(i)+0.0098*blk_ZT(ng))
# endif
# ifdef ICE_DIAGS
sens_ice(i,j) = qsh_i
# endif
! Latent heat of sublimation
le_i = 2.834E+6_r8
# ifdef ICE_BOX
qlh_i = latent_MU(month)
# ifdef ICE_I_O
SW_thru_ice(i,j) = 0.0_r8
# endif
# else
! Qsati is saturation specific humidity over the ice from
! Parkinson and Washington (1979)
slp = Pair(i,j)*100._r8 !Convert back to Pascal from millibars
cff = 611._r8* &
& 10._r8**(9.5_r8*(sfc_temp-273.16_r8)/(sfc_temp-7.66_r8))
Qsati=epsilon*cff/(slp-(1._r8-epsilon)*cff)
dq_i=Hair(i,j)-Qsati
fqlat1=rhoAir(i)*Ce*vmag
qlh_i = fqlat1*le_i*dq_i
# endif
# ifdef ICE_DIAGS
lat_ice(i,j) = qlh_i
# endif