00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 module ice_forcing
00028
00029
00030
00031 use ice_kinds_mod
00032 use ice_blocks, only: nx_block, ny_block
00033 use ice_domain_size
00034 use ice_communicate, only: my_task, master_task
00035 use ice_constants
00036 use ice_calendar, only: istep, istep1, time, time_forc, year_init, &
00037 sec, mday, month, nyr, yday, daycal, dayyr, &
00038 daymo, days_per_year
00039 use ice_fileunits
00040 use ice_atmo, only: calc_strair
00041 use ice_exit
00042 use ice_timers
00043
00044
00045
00046 implicit none
00047 save
00048
00049 integer (kind=int_kind) ::
00050 ycycle ,
00051 fyear_init ,
00052 fyear ,
00053 fyear_final
00054
00055 character (char_len_long) ::
00056 height_file,
00057 uwind_file,
00058 vwind_file,
00059 wind_file,
00060 strax_file,
00061 stray_file,
00062 potT_file,
00063 tair_file,
00064 humid_file,
00065 rhoa_file,
00066 fsw_file,
00067 flw_file,
00068 rain_file,
00069 sst_file,
00070 sss_file,
00071 pslv_file,
00072 sublim_file,
00073 snow_file
00074
00075 character (char_len_long), dimension(ncat) ::
00076 topmelt_file,
00077 botmelt_file
00078
00079 real (kind=dbl_kind) ::
00080 c1intp, c2intp ,
00081 ftime
00082
00083 integer (kind=int_kind) ::
00084 oldrecnum = 0 ,
00085 oldrecslot = 1
00086
00087 real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) ::
00088 cldf
00089
00090 real (kind=dbl_kind), dimension(nx_block,ny_block,2,max_blocks) ::
00091 fsw_data,
00092 cldf_data,
00093 fsnow_data,
00094 Tair_data,
00095 uatm_data,
00096 vatm_data,
00097 wind_data,
00098 strax_data,
00099 stray_data,
00100 Qa_data,
00101 rhoa_data,
00102 potT_data,
00103 zlvl_data,
00104 flw_data,
00105 sst_data,
00106 sss_data,
00107 uocn_data,
00108 vocn_data,
00109 sublim_data,
00110 frain_data
00111
00112 real (kind=dbl_kind),
00113 dimension(nx_block,ny_block,2,max_blocks,ncat) ::
00114 topmelt_data,
00115 botmelt_data
00116
00117 character(char_len) ::
00118 atm_data_format,
00119 ocn_data_format,
00120 atm_data_type,
00121
00122 sss_data_type, &
00123 sst_data_type, &
00124
00125 precip_units
00126
00127 character(char_len_long) ::
00128 atm_data_dir ,
00129 ocn_data_dir ,
00130 oceanmixed_file
00131
00132 integer (kind=int_kind), parameter ::
00133 nfld = 8
00134
00135
00136 real (kind=dbl_kind), parameter ::
00137 frcvdr = 0.28_dbl_kind,
00138 frcvdf = 0.24_dbl_kind,
00139 frcidr = 0.31_dbl_kind,
00140 frcidf = 0.17_dbl_kind
00141
00142 real (kind=dbl_kind),
00143 dimension (nx_block,ny_block,max_blocks,nfld,12) ::
00144 ocn_frc_m
00145
00146 logical (kind=log_kind) ::
00147 restore_sst
00148
00149 integer (kind=int_kind) ::
00150 trestore
00151
00152 real (kind=dbl_kind) ::
00153 trest
00154
00155 logical (kind=log_kind) ::
00156 dbug
00157
00158
00159
00160 contains
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170 subroutine init_forcing_atmo
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187 fyear = fyear_init + mod(nyr-1,ycycle)
00188 fyear_final = fyear_init + ycycle - 1
00189
00190 if (trim(atm_data_type) /= 'default' .and. &
00191 my_task == master_task) then
00192 write (nu_diag,*) ' Initial forcing data year = ',fyear_init
00193 write (nu_diag,*) ' Final forcing data year = ',fyear_final
00194 endif
00195
00196
00197
00198
00199
00200
00201 if (trim(atm_data_type) == 'ncar') then
00202 call NCAR_files(fyear)
00203 elseif (trim(atm_data_type) == 'ecmwf') then
00204 call ecmwf_files(fyear)
00205 elseif (trim(atm_data_type) == 'LYq') then
00206 call LY_files(fyear)
00207 elseif (trim(atm_data_type) == 'hadgem') then
00208 call hadgem_files(fyear)
00209 elseif (trim(atm_data_type) == 'monthly') then
00210 call monthly_files(fyear)
00211 endif
00212
00213 end subroutine init_forcing_atmo
00214
00215
00216
00217
00218
00219
00220
00221
00222 subroutine init_forcing_ocn(dt)
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241 use ice_domain, only: nblocks
00242 use ice_flux, only: sss, sst, Tf, Tfrzpt
00243 use ice_work, only: work1
00244 use ice_read_write
00245 use ice_therm_vertical, only: ustar_scale
00246
00247
00248
00249 real (kind=dbl_kind), intent(in) ::
00250 dt
00251
00252
00253
00254 integer (kind=int_kind) ::
00255 i, j, iblk ,
00256 k ,
00257 fid ,
00258 nbits
00259
00260 logical (kind=log_kind) :: diag
00261
00262 character (char_len) ::
00263 fieldname
00264
00265 nbits = 64
00266
00267
00268
00269
00270
00271
00272 if (trim(sss_data_type) == 'clim') then
00273
00274
00275 sss_file = trim(ocn_data_dir)//'sss.mm.100x116.da'
00276
00277
00278 if (my_task == master_task) then
00279 write (nu_diag,*) ' '
00280 write (nu_diag,*) 'SSS climatology computed from:'
00281 write (nu_diag,*) trim(sss_file)
00282 endif
00283
00284 if (my_task == master_task) &
00285 call ice_open (nu_forcing, sss_file, nbits)
00286
00287 sss(:,:,:) = c0
00288
00289 do k = 1,12
00290 call ice_read (nu_forcing, k, work1, 'rda8', dbug, &
00291 field_loc_center, field_type_scalar)
00292 do iblk = 1, nblocks
00293 do j = 1, ny_block
00294 do i = 1, nx_block
00295 sss(i,j,iblk) = sss(i,j,iblk) + work1(i,j,iblk)
00296 enddo
00297 enddo
00298 enddo
00299 enddo
00300
00301 do iblk = 1, nblocks
00302 do j = 1, ny_block
00303 do i = 1, nx_block
00304 sss(i,j,iblk) = sss(i,j,iblk) / c12
00305 sss(i,j,iblk) = max(sss(i,j,iblk),c0)
00306 if (trim(Tfrzpt) == 'constant') then
00307 Tf (i,j,iblk) = -1.8_dbl_kind
00308 else
00309 Tf (i,j,iblk) = -depressT * sss(i,j,iblk)
00310 endif
00311 enddo
00312 enddo
00313 enddo
00314
00315
00316 if (my_task == master_task) close(nu_forcing)
00317
00318 endif
00319
00320
00321
00322
00323
00324
00325 if (restore_sst) then
00326 if (trestore == 0) then
00327 trest = dt
00328 else
00329 trest = real(trestore,kind=dbl_kind) * secday
00330 endif
00331 endif
00332
00333 if (trim(sst_data_type) == 'clim') then
00334
00335 if (nx_global == 320) then
00336 sst_file = trim(ocn_data_dir)//'sst_clim_hurrell.dat'
00337 else
00338
00339 sst_file = trim(ocn_data_dir)//'sst.mm.100x116.da'
00340
00341 endif
00342
00343 if (my_task == master_task) then
00344 write (nu_diag,*) ' '
00345 write (nu_diag,*) 'Initial SST file:', trim(sst_file)
00346 endif
00347
00348 if (my_task == master_task) &
00349 call ice_open (nu_forcing, sst_file, nbits)
00350
00351 call ice_read (nu_forcing, month, sst, 'rda8', dbug, &
00352 field_loc_center, field_type_scalar)
00353
00354 if (my_task == master_task) close(nu_forcing)
00355
00356
00357 do iblk = 1, nblocks
00358 do j = 1, ny_block
00359 do i = 1, nx_block
00360 sst(i,j,iblk) = max(sst(i,j,iblk),Tf(i,j,iblk))
00361 enddo
00362 enddo
00363 enddo
00364
00365 endif
00366
00367
00368 if (trim(sst_data_type) == 'hadgem_sst' .or. &
00369 trim(sst_data_type) == 'hadgem_sst_uvocn') then
00370
00371 diag = .true.
00372
00373 sst_file = trim (ocn_data_dir)//'MONTHLY/sst.1997.nc'
00374
00375 if (my_task == master_task) then
00376
00377 write (nu_diag,*) ' '
00378 write (nu_diag,*) 'Initial SST file:', trim(sst_file)
00379
00380 call ice_open_nc(sst_file,fid)
00381
00382 endif
00383
00384 fieldname='sst'
00385 call ice_read_nc(fid,month,fieldname,sst,diag)
00386
00387 if (my_task == master_task) call ice_close_nc(fid)
00388
00389
00390 do iblk = 1, nblocks
00391 do j = 1, ny_block
00392 do i = 1, nx_block
00393 sst(i,j,iblk) = max(sst(i,j,iblk),Tf(i,j,iblk))
00394 enddo
00395 enddo
00396 enddo
00397
00398 endif
00399
00400
00401 if (trim(sst_data_type) == 'ncar' .or. &
00402 trim(sss_data_type) == 'ncar') then
00403 call ocn_data_ncar_init
00404 endif
00405
00406
00407
00408
00409 if (trim(sst_data_type) /= 'ncar' .or. &
00410 trim(sss_data_type) /= 'ncar' .or. &
00411 trim(sst_data_type) /= 'hadgem_sst_uvocn') then
00412 ustar_scale = c10
00413 endif
00414
00415 end subroutine init_forcing_ocn
00416
00417
00418
00419
00420
00421
00422
00423
00424 subroutine get_forcing_atmo
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436 use ice_boundary, only: ice_HaloUpdate
00437 use ice_domain
00438 use ice_blocks
00439 use ice_flux
00440 use ice_state
00441 use ice_grid, only: ANGLET, hm
00442
00443
00444
00445 integer (kind=int_kind) ::
00446 iblk,
00447 ilo,ihi,jlo,jhi
00448
00449 type (block) ::
00450 this_block
00451
00452 fyear = fyear_init + mod(nyr-1,ycycle)
00453 if (trim(atm_data_type) /= 'default' .and. istep <= 1 &
00454 .and. my_task == master_task) then
00455 write (nu_diag,*) ' Current forcing data year = ',fyear
00456 endif
00457
00458 ftime = time
00459 time_forc = ftime
00460
00461
00462
00463
00464
00465 if (trim(atm_data_type) == 'ncar') then
00466 call ncar_data
00467 elseif (trim(atm_data_type) == 'ecmwf') then
00468 call ecmwf_data
00469 elseif (trim(atm_data_type) == 'LYq') then
00470 call LY_data
00471 elseif (trim(atm_data_type) == 'hadgem') then
00472 call hadgem_data
00473 elseif (trim(atm_data_type) == 'monthly') then
00474 call monthly_data
00475 else
00476 return
00477 endif
00478
00479
00480
00481
00482
00483 do iblk = 1, nblocks
00484
00485 this_block = get_block(blocks_ice(iblk),iblk)
00486 ilo = this_block%ilo
00487 ihi = this_block%ihi
00488 jlo = this_block%jlo
00489 jhi = this_block%jhi
00490
00491 call prepare_forcing (nx_block, ny_block, &
00492 ilo, ihi, jlo, jhi, &
00493 hm (:,:,iblk), &
00494 Tair (:,:,iblk), &
00495 fsw (:,:,iblk), &
00496 cldf (:,:,iblk), &
00497 flw (:,:,iblk), &
00498 frain (:,:,iblk), &
00499 fsnow (:,:,iblk), &
00500 Qa (:,:,iblk), &
00501 rhoa (:,:,iblk), &
00502 uatm (:,:,iblk), &
00503 vatm (:,:,iblk), &
00504 strax (:,:,iblk), &
00505 stray (:,:,iblk), &
00506 zlvl (:,:,iblk), &
00507 wind (:,:,iblk), &
00508 swvdr (:,:,iblk), &
00509 swvdf (:,:,iblk), &
00510 swidr (:,:,iblk), &
00511 swidf (:,:,iblk), &
00512 potT (:,:,iblk), &
00513 ANGLET(:,:,iblk), &
00514 trcr (:,:,nt_Tsfc,iblk), &
00515 sst (:,:,iblk), &
00516 aice (:,:,iblk) )
00517
00518 enddo
00519
00520 call ice_timer_start(timer_bound)
00521 call ice_HaloUpdate (swvdr, halo_info, &
00522 field_loc_center, field_type_scalar)
00523 call ice_HaloUpdate (swvdf, halo_info, &
00524 field_loc_center, field_type_scalar)
00525 call ice_HaloUpdate (swidr, halo_info, &
00526 field_loc_center, field_type_scalar)
00527 call ice_HaloUpdate (swidf, halo_info, &
00528 field_loc_center, field_type_scalar)
00529 call ice_timer_stop(timer_bound)
00530
00531 end subroutine get_forcing_atmo
00532
00533
00534
00535
00536
00537
00538
00539
00540 subroutine get_forcing_ocn (dt)
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554 use ice_domain, only: nblocks
00555 use ice_ocean
00556
00557
00558
00559 real (kind=dbl_kind), intent(in) ::
00560 dt
00561
00562
00563
00564 if (trim(sst_data_type) == 'clim' .or. &
00565 trim(sss_data_type) == 'clim') then
00566 call ocn_data_clim(dt)
00567 elseif (trim(sst_data_type) == 'ncar' .or. &
00568 trim(sss_data_type) == 'ncar') then
00569 call ocn_data_ncar(dt)
00570 elseif (trim(sst_data_type) == 'hadgem_sst' .or. &
00571 trim(sst_data_type) == 'hadgem_sst_uvocn') then
00572 call ocn_data_hadgem(dt)
00573 endif
00574
00575 end subroutine get_forcing_ocn
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585 subroutine read_data (flag, recd, yr, ixm, ixx, ixp, &
00586 maxrec, data_file, field_data, &
00587 field_loc, field_type)
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613 use ice_read_write
00614 use ice_diagnostics, only: check_step
00615
00616
00617
00618 logical (kind=log_kind), intent(in) :: flag
00619
00620 integer (kind=int_kind), intent(in) ::
00621 recd ,
00622 yr ,
00623 ixm, ixx, ixp ,
00624
00625 maxrec
00626
00627 real (kind=dbl_kind), dimension(nx_block,ny_block,2,max_blocks),
00628 intent(out) ::
00629 field_data
00630
00631 integer (kind=int_kind), intent(in) ::
00632 field_loc,
00633 field_type
00634
00635
00636
00637 character (char_len_long) ::
00638 data_file
00639
00640 integer (kind=int_kind) ::
00641 nbits ,
00642 nrec ,
00643 n2, n4 ,
00644
00645 arg
00646
00647 call ice_timer_start(timer_readwrite)
00648
00649 nbits = 64
00650
00651 if (istep1 > check_step) dbug = .true.
00652
00653 if (my_task==master_task .and. (dbug)) then
00654 write(nu_diag,*) ' ', trim(data_file)
00655 endif
00656
00657 if (flag) then
00658
00659
00660
00661
00662
00663
00664 n2 = ixm
00665 n4 = ixp
00666 arg = 0
00667
00668
00669
00670
00671
00672 if (ixm /= 99) then
00673
00674 if (ixx <= 1) then
00675 if (yr > fyear_init) then
00676 call file_year (data_file, yr-1)
00677 else
00678 if (maxrec > 12) then
00679 if (ixx == 1) n2 = ixx
00680 else
00681 call file_year (data_file, fyear_final)
00682 endif
00683 endif
00684 endif
00685
00686 call ice_open (nu_forcing, data_file, nbits)
00687
00688 arg = 1
00689 nrec = recd + n2
00690 call ice_read (nu_forcing, nrec, field_data(:,:,arg,:), &
00691 'rda8', dbug, field_loc, field_type)
00692
00693 if (ixx==1 .and. my_task == master_task) close(nu_forcing)
00694 endif
00695
00696
00697 call file_year (data_file, yr)
00698 call ice_open (nu_forcing, data_file, nbits)
00699
00700 arg = arg + 1
00701 nrec = recd + ixx
00702 call ice_read (nu_forcing, nrec, field_data(:,:,arg,:), &
00703 'rda8', dbug, field_loc, field_type)
00704
00705 if (ixp /= 99) then
00706
00707 if (ixx==maxrec) then
00708 if (yr < fyear_final) then
00709 if (my_task == master_task) close(nu_forcing)
00710 call file_year (data_file, yr+1)
00711 call ice_open (nu_forcing, data_file, nbits)
00712 else
00713 if (maxrec > 12) then
00714 n4 = ixx
00715 else
00716 if (my_task == master_task) close(nu_forcing)
00717 call file_year (data_file, fyear_init)
00718
00719 call ice_open (nu_forcing, data_file, nbits)
00720
00721 endif
00722 endif
00723 endif
00724
00725 arg = arg + 1
00726 nrec = recd + n4
00727 call ice_read (nu_forcing, nrec, field_data(:,:,arg,:), &
00728 'rda8', dbug, field_loc, field_type)
00729 endif
00730
00731 if (my_task == master_task) close(nu_forcing)
00732
00733 endif
00734
00735 call ice_timer_stop(timer_readwrite)
00736
00737 end subroutine read_data
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747 subroutine read_data_nc (flag, recd, yr, ixm, ixx, ixp, &
00748 maxrec, data_file, fieldname, field_data, &
00749 field_loc, field_type)
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775 use ice_read_write
00776 use ice_diagnostics, only: check_step
00777
00778
00779
00780 logical (kind=log_kind), intent(in) :: flag
00781
00782 integer (kind=int_kind), intent(in) ::
00783 recd ,
00784 yr ,
00785 ixm, ixx, ixp ,
00786
00787 maxrec
00788
00789 character (char_len_long) ::
00790 data_file
00791
00792 character (char_len), intent(in) ::
00793 fieldname
00794
00795 integer (kind=int_kind), intent(in) ::
00796 field_loc,
00797 field_type
00798
00799 real (kind=dbl_kind), dimension(nx_block,ny_block,2,max_blocks),
00800 intent(out) ::
00801 field_data
00802
00803
00804
00805 #ifdef ncdf
00806 integer (kind=int_kind) ::
00807 nrec ,
00808 n2, n4 ,
00809
00810 arg , &
00811 fid
00812
00813
00814 call ice_timer_start(timer_readwrite)
00815
00816 if (istep1 > check_step) dbug = .true.
00817
00818 if (my_task==master_task .and. (dbug)) then
00819 write(nu_diag,*) ' ', trim(data_file)
00820 endif
00821
00822 if (flag) then
00823
00824
00825
00826
00827
00828
00829 n2 = ixm
00830 n4 = ixp
00831 arg = 0
00832
00833
00834
00835
00836
00837 if (ixm /= 99) then
00838
00839 if (ixx <= 1) then
00840 if (yr > fyear_init) then
00841 call file_year (data_file, yr-1)
00842 else
00843 if (maxrec > 12) then
00844 if (ixx == 1) n2 = ixx
00845 else
00846 call file_year (data_file, fyear_final)
00847 endif
00848 endif
00849 endif
00850
00851 call ice_open_nc (data_file, fid)
00852
00853 arg = 1
00854 nrec = recd + n2
00855
00856 call ice_read_nc &
00857 (fid, nrec, fieldname, field_data(:,:,arg,:), dbug, &
00858 field_loc, field_type)
00859
00860 if (ixx==1) call ice_close_nc(fid)
00861 endif
00862
00863
00864 call file_year (data_file, yr)
00865 call ice_open_nc (data_file, fid)
00866
00867 arg = arg + 1
00868 nrec = recd + ixx
00869
00870 call ice_read_nc &
00871 (fid, nrec, fieldname, field_data(:,:,arg,:), dbug, &
00872 field_loc, field_type)
00873
00874 if (ixp /= 99) then
00875
00876 if (ixx==maxrec) then
00877 if (yr < fyear_final) then
00878 call ice_close_nc(fid)
00879 call file_year (data_file, yr+1)
00880 call ice_open_nc (data_file, fid)
00881 else
00882 if (maxrec > 12) then
00883 n4 = ixx
00884 else
00885 call ice_close_nc(fid)
00886 call file_year (data_file, fyear_init)
00887 call ice_open_nc (data_file, fid)
00888
00889 endif
00890 endif
00891 endif
00892
00893 arg = arg + 1
00894 nrec = recd + n4
00895
00896 call ice_read_nc &
00897 (fid, nrec, fieldname, field_data(:,:,arg,:), dbug, &
00898 field_loc, field_type)
00899 endif
00900
00901 call ice_close_nc(fid)
00902
00903 endif
00904
00905 call ice_timer_stop(timer_readwrite)
00906
00907 #else
00908 field_data = c0
00909 #endif
00910 end subroutine read_data_nc
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920 subroutine read_clim_data (readflag, recd, ixm, ixx, ixp, &
00921 data_file, field_data, &
00922 field_loc, field_type)
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937 use ice_read_write
00938 use ice_diagnostics, only: check_step
00939
00940
00941
00942 logical (kind=log_kind),intent(in) :: readflag
00943
00944 integer (kind=int_kind), intent(in) ::
00945 recd ,
00946 ixm,ixx,ixp
00947
00948
00949 character (char_len_long), intent(in) :: data_file
00950
00951 integer (kind=int_kind), intent(in) ::
00952 field_loc,
00953 field_type
00954
00955 real (kind=dbl_kind), dimension(nx_block,ny_block,2,max_blocks),
00956 intent(out) ::
00957 field_data
00958
00959
00960
00961 integer (kind=int_kind) ::
00962 nbits ,
00963 nrec ,
00964 arg
00965
00966 call ice_timer_start(timer_readwrite)
00967
00968 nbits = 64
00969
00970 if (istep1 > check_step) dbug = .true.
00971
00972 if (my_task==master_task .and. (dbug)) &
00973 write(nu_diag,*) ' ', trim(data_file)
00974
00975 if (readflag) then
00976
00977
00978
00979
00980
00981 call ice_open (nu_forcing, data_file, nbits)
00982
00983 arg = 0
00984 if (ixm /= 99) then
00985 arg = 1
00986 nrec = recd + ixm
00987 call ice_read (nu_forcing, nrec, field_data(:,:,arg,:), &
00988 'rda8', dbug, field_loc, field_type)
00989 endif
00990
00991 arg = arg + 1
00992 nrec = recd + ixx
00993 call ice_read (nu_forcing, nrec, field_data(:,:,arg,:), &
00994 'rda8', dbug, field_loc, field_type)
00995
00996 if (ixp /= 99) then
00997 arg = arg + 1
00998 nrec = recd + ixp
00999 call ice_read (nu_forcing, nrec, field_data(:,:,arg,:), &
01000 'rda8', dbug, field_loc, field_type)
01001 endif
01002
01003 if (my_task == master_task) close (nu_forcing)
01004 endif
01005
01006 call ice_timer_stop(timer_readwrite)
01007
01008 end subroutine read_clim_data
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018 subroutine interp_coeff_monthly (recslot)
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032 integer (kind=int_kind), intent(in) ::
01033 recslot
01034
01035
01036
01037 real (kind=dbl_kind) ::
01038 tt ,
01039 t1, t2
01040
01041 real (kind=dbl_kind) ::
01042 daymid(0:13)
01043
01044 daymid(1:13) = 14._dbl_kind
01045 daymid(0) = 14._dbl_kind - daymo(12)
01046
01047
01048 tt = mod(ftime/secday,dayyr)
01049
01050
01051
01052 if (recslot==2) then
01053 t2 = daycal(month) + daymid(month)
01054 if (month == 1) then
01055 t1 = daymid(0)
01056 else
01057 t1 = daycal(month-1) + daymid(month-1)
01058 endif
01059 else
01060 t1 = daycal(month) + daymid(month)
01061 t2 = daycal(month+1) + daymid(month+1)
01062 endif
01063
01064
01065 c1intp = (t2 - tt) / (t2 - t1)
01066 c2intp = c1 - c1intp
01067
01068 end subroutine interp_coeff_monthly
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078 subroutine interp_coeff (recnum, recslot, secint, dataloc)
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093 integer (kind=int_kind), intent(in) ::
01094 recnum ,
01095 recslot ,
01096 dataloc
01097
01098
01099 real (kind=dbl_kind), intent(in) ::
01100 secint
01101
01102
01103
01104
01105
01106 real (kind=dbl_kind) ::
01107 secyr
01108
01109 real (kind=dbl_kind) ::
01110 tt ,
01111 t1, t2 ,
01112 rcnum
01113
01114 secyr = dayyr * secday
01115 tt = mod(ftime,secyr)
01116
01117
01118 rcnum = real(recnum,kind=dbl_kind)
01119 if (recslot==2) then
01120 if (dataloc==1) then
01121 t2 = (rcnum-p5)*secint
01122 else
01123 t2 = rcnum*secint
01124 endif
01125 t1 = t2 - secint
01126 else
01127 if (dataloc==1) then
01128 t1 = (rcnum-p5)*secint
01129 else
01130 t1 = rcnum*secint
01131 endif
01132 t2 = t1 + secint
01133 endif
01134
01135
01136 c1intp = abs((t2 - tt) / (t2 - t1))
01137 c2intp = c1 - c1intp
01138
01139 end subroutine interp_coeff
01140
01141
01142
01143
01144
01145
01146
01147
01148 subroutine interpolate_data (field_data, field)
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160 use ice_domain, only: nblocks
01161
01162
01163
01164 real (kind=dbl_kind), dimension(nx_block,ny_block,2,max_blocks),
01165 intent(in) ::
01166 field_data
01167
01168 real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks),
01169 intent(out) ::
01170 field
01171
01172
01173
01174 integer (kind=int_kind) :: i,j, iblk
01175
01176 do iblk = 1, nblocks
01177 do j = 1, ny_block
01178 do i = 1, nx_block
01179 field(i,j,iblk) = c1intp * field_data(i,j,1,iblk) &
01180 + c2intp * field_data(i,j,2,iblk)
01181 enddo
01182 enddo
01183 enddo
01184
01185 end subroutine interpolate_data
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195 subroutine file_year (data_file, yr)
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212 character (char_len_long), intent(inout) :: data_file
01213
01214
01215
01216 integer (kind=int_kind), intent(in) :: yr
01217
01218 character (char_len_long) :: tmpname
01219
01220 integer (kind=int_kind) :: i
01221
01222 if (trim(atm_data_type) == 'ecmwf') then
01223 i = index(data_file,'.r') - 5
01224 tmpname = data_file
01225 write(data_file,'(a,i4.4,a)') tmpname(1:i), yr, '.r'
01226 elseif (trim(atm_data_type) == 'hadgem') then
01227 i = index(data_file,'.nc') - 5
01228 tmpname = data_file
01229 write(data_file,'(a,i4.4,a)') tmpname(1:i), yr, '.nc'
01230 else
01231 i = index(data_file,'.dat') - 5
01232 tmpname = data_file
01233 write(data_file,'(a,i4.4,a)') tmpname(1:i), yr, '.dat'
01234 endif
01235
01236 end subroutine file_year
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246 subroutine prepare_forcing (nx_block, ny_block, &
01247 ilo, ihi, jlo, jhi, &
01248 hm, &
01249 Tair, fsw, &
01250 cldf, flw, &
01251 frain, fsnow, &
01252 Qa, rhoa, &
01253 uatm, vatm, &
01254 strax, stray, &
01255 zlvl, wind, &
01256 swvdr, swvdf, &
01257 swidr, swidf, &
01258 potT, ANGLET, &
01259 Tsfc, sst, &
01260 aice)
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272 integer (kind=int_kind), intent(in) ::
01273 nx_block, ny_block,
01274 ilo,ihi,jlo,jhi
01275
01276 real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) ::
01277 Tair ,
01278 ANGLET ,
01279 Tsfc ,
01280 sst ,
01281 aice ,
01282 hm
01283
01284 real (kind=dbl_kind), dimension(nx_block,ny_block),
01285 intent(inout) ::
01286 fsw ,
01287 cldf ,
01288 frain ,
01289 fsnow ,
01290 Qa ,
01291 rhoa ,
01292 uatm ,
01293 vatm ,
01294 strax ,
01295 stray ,
01296 zlvl ,
01297 wind ,
01298 flw ,
01299 swvdr ,
01300 swvdf ,
01301 swidr ,
01302 swidf ,
01303 potT
01304
01305
01306
01307 integer (kind=int_kind) ::
01308 i, j
01309
01310 real (kind=dbl_kind) :: workx, worky,
01311 fcc, sstk, rtea, ptem, qlwm, precip_factor
01312
01313 do j = jlo, jhi
01314 do i = ilo, ihi
01315
01316
01317
01318
01319 cldf (i,j) = max(min(cldf(i,j),c1),c0)
01320 fsw (i,j) = max(fsw(i,j),c0)
01321 fsnow(i,j) = max(fsnow(i,j),c0)
01322 rhoa (i,j) = max(rhoa(i,j),c0)
01323 Qa (i,j) = max(Qa(i,j),c0)
01324
01325 enddo
01326 enddo
01327
01328
01329
01330
01331
01332 if (trim(atm_data_type) == 'ncar') then
01333 do j = jlo, jhi
01334 do i = ilo, ihi
01335
01336
01337
01338
01339
01340 Qa (i,j) = Qa (i,j) * 0.94_dbl_kind
01341 fsw(i,j) = fsw(i,j) * 0.92_dbl_kind
01342
01343
01344
01345
01346
01347
01348 flw(i,j) = stefan_boltzmann*Tair(i,j)**4 &
01349 * (c1 - 0.261_dbl_kind &
01350 *exp(-7.77e-4_dbl_kind*(Tffresh - Tair(i,j))**2)) &
01351 * (c1 + 0.275_dbl_kind*cldf(i,j))
01352
01353
01354
01355 enddo
01356 enddo
01357
01358 elseif (trim(atm_data_type) == 'ecmwf') then
01359 do j = jlo, jhi
01360 do i = ilo, ihi
01361
01362
01363
01364
01365
01366
01367
01368 Qa (i,j) = (qqqocn/rhoa(i,j)) * exp(-TTTocn/Qa(i,j))
01369 Qa (i,j) = max(Qa(i,j),c0)
01370
01371
01372
01373
01374
01375
01376 flw(i,j) = stefan_boltzmann*Tair(i,j)**4 &
01377 * (c1 - 0.261_dbl_kind &
01378 *exp(-7.77e-4_dbl_kind*(Tffresh - Tair(i,j))**2)) &
01379 * (c1 + 0.275_dbl_kind*cldf(i,j))
01380
01381
01382
01383 enddo
01384 enddo
01385
01386 elseif (trim(atm_data_type) == 'LYq') then
01387
01388
01389
01390
01391
01392 do j = jlo, jhi
01393 do i = ilo, ihi
01394 fcc = c1 - 0.8_dbl_kind * cldf(i,j)
01395 sstk = (Tsfc(i,j) * aice(i,j) &
01396 + sst(i,j) * (c1 - aice(i,j))) + Tffresh
01397 rtea = sqrt(c1000*Qa(i,j) / &
01398 (0.622_dbl_kind+0.378_dbl_kind*Qa(i,j)))
01399 ptem = Tair(i,j)
01400 qlwm = ptem * ptem * ptem &
01401 * ( ptem*(0.39_dbl_kind-0.05_dbl_kind*rtea)*fcc &
01402 + c4*(sstk-ptem) )
01403 flw(i,j) = emissivity*stefan_boltzmann * ( sstk**4 - qlwm )
01404 flw(i,j) = flw(i,j) * hm(i,j)
01405 enddo
01406 enddo
01407
01408 endif
01409
01410
01411
01412
01413
01414
01415 if (trim(precip_units) == 'mm_per_month') then
01416 precip_factor = c12/(secday*days_per_year)
01417 elseif (trim(precip_units) == 'mm_per_day') then
01418 precip_factor = c1/secday
01419 elseif (trim(precip_units) == 'mm_per_sec' .or. &
01420 trim(precip_units) == 'mks') then
01421 precip_factor = c1
01422 endif
01423
01424 do j = jlo, jhi
01425 do i = ilo, ihi
01426
01427 zlvl(i,j) = c10
01428 potT(i,j) = Tair(i,j)
01429
01430
01431 swvdr(i,j) = fsw(i,j)*frcvdr
01432 swvdf(i,j) = fsw(i,j)*frcvdf
01433 swidr(i,j) = fsw(i,j)*frcidr
01434 swidf(i,j) = fsw(i,j)*frcidf
01435
01436
01437 fsnow(i,j) = fsnow(i,j) * precip_factor
01438
01439 enddo
01440 enddo
01441
01442
01443
01444
01445 if (trim(atm_data_type) /= 'hadgem') then
01446
01447 do j = jlo, jhi
01448 do i = ilo, ihi
01449 frain(i,j) = c0
01450 if (Tair(i,j) >= Tffresh) then
01451 frain(i,j) = fsnow(i,j)
01452 fsnow(i,j) = c0
01453 endif
01454 enddo
01455 enddo
01456
01457 endif
01458
01459 if (calc_strair) then
01460
01461 do j = jlo, jhi
01462 do i = ilo, ihi
01463
01464 wind(i,j) = sqrt(uatm(i,j)**2 + vatm(i,j)**2)
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477 workx = uatm(i,j)
01478 worky = vatm(i,j)
01479 uatm (i,j) = workx*cos(ANGLET(i,j)) &
01480 + worky*sin(ANGLET(i,j))
01481 vatm (i,j) = worky*cos(ANGLET(i,j)) &
01482 - workx*sin(ANGLET(i,j))
01483
01484 enddo
01485 enddo
01486
01487 else
01488
01489 do j = jlo, jhi
01490 do i = ilo, ihi
01491
01492 workx = strax(i,j)
01493 worky = stray(i,j)
01494 strax(i,j) = workx*cos(ANGLET(i,j)) &
01495 + worky*sin(ANGLET(i,j))
01496 stray(i,j) = worky*cos(ANGLET(i,j)) &
01497 - workx*sin(ANGLET(i,j))
01498
01499 enddo
01500 enddo
01501
01502 endif
01503
01504 end subroutine prepare_forcing
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516 subroutine ncar_files (yr)
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534 integer (kind=int_kind), intent(in) ::
01535 yr
01536
01537
01538
01539 fsw_file = &
01540 trim(atm_data_dir)//'ISCCPM/MONTHLY/RADFLX/swdn.1996.dat'
01541 call file_year(fsw_file,yr)
01542
01543 flw_file = &
01544 trim(atm_data_dir)//'ISCCPM/MONTHLY/RADFLX/cldf.1996.dat'
01545 call file_year(flw_file,yr)
01546
01547 rain_file = &
01548 trim(atm_data_dir)//'MXA/MONTHLY/PRECIP/prec.1996.dat'
01549 call file_year(rain_file,yr)
01550
01551 uwind_file = &
01552 trim(atm_data_dir)//'NCEP/4XDAILY/STATES/u_10.1996.dat'
01553 call file_year(uwind_file,yr)
01554
01555 vwind_file = &
01556 trim(atm_data_dir)//'NCEP/4XDAILY/STATES/v_10.1996.dat'
01557 call file_year(vwind_file,yr)
01558
01559 tair_file = &
01560 trim(atm_data_dir)//'NCEP/4XDAILY/STATES/t_10.1996.dat'
01561 call file_year(tair_file,yr)
01562
01563 humid_file = &
01564 trim(atm_data_dir)//'NCEP/4XDAILY/STATES/q_10.1996.dat'
01565 call file_year(humid_file,yr)
01566
01567 rhoa_file = &
01568 trim(atm_data_dir)//'NCEP/4XDAILY/STATES/dn10.1996.dat'
01569 call file_year(rhoa_file,yr)
01570
01571 if (my_task == master_task) then
01572 write (nu_diag,*) ' '
01573 write (nu_diag,*) 'Forcing data year =', fyear
01574 write (nu_diag,*) 'Atmospheric data files:'
01575 write (nu_diag,*) trim(fsw_file)
01576 write (nu_diag,*) trim(flw_file)
01577 write (nu_diag,*) trim(rain_file)
01578 write (nu_diag,*) trim(uwind_file)
01579 write (nu_diag,*) trim(vwind_file)
01580 write (nu_diag,*) trim(tair_file)
01581 write (nu_diag,*) trim(humid_file)
01582 write (nu_diag,*) trim(rhoa_file)
01583 endif
01584
01585 end subroutine ncar_files
01586
01587
01588
01589
01590
01591
01592
01593
01594
01595 subroutine ncar_data
01596
01597
01598
01599
01600
01601
01602
01603
01604
01605 use ice_flux
01606
01607
01608
01609
01610
01611
01612 integer (kind=int_kind) ::
01613 i, j ,
01614 ixm,ixx,ixp ,
01615 recnum ,
01616 maxrec ,
01617 recslot ,
01618 dataloc ,
01619
01620 midmonth
01621
01622 real (kind=dbl_kind) ::
01623 sec6hr
01624
01625 logical (kind=log_kind) :: readm, read6
01626
01627
01628
01629
01630
01631
01632
01633
01634 midmonth = 15
01635
01636
01637
01638 maxrec = 12
01639 ixm = mod(month+maxrec-2,maxrec) + 1
01640 ixp = mod(month, maxrec) + 1
01641 if (mday >= midmonth) ixm = 99
01642 if (mday < midmonth) ixp = 99
01643
01644
01645
01646
01647
01648
01649 recslot = 1
01650 if (mday < midmonth) recslot = 2
01651
01652
01653 call interp_coeff_monthly (recslot)
01654
01655
01656 readm = .false.
01657 if (istep==1 .or. (mday==midmonth .and. sec==0)) readm = .true.
01658
01659 if (trim(atm_data_format) == 'bin') then
01660 call read_data (readm, 0, fyear, ixm, month, ixp, &
01661 maxrec, fsw_file, fsw_data, &
01662 field_loc_center, field_type_scalar)
01663 call read_data (readm, 0, fyear, ixm, month, ixp, &
01664 maxrec, flw_file, cldf_data, &
01665 field_loc_center, field_type_scalar)
01666 call read_data (readm, 0, fyear, ixm, month, ixp, &
01667 maxrec, rain_file, fsnow_data, &
01668 field_loc_center, field_type_scalar)
01669 else
01670 call abort_ice ('nonbinary atm_data_format unavailable')
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681 endif
01682
01683
01684 call interpolate_data (fsw_data, fsw)
01685 call interpolate_data (cldf_data, cldf)
01686 call interpolate_data (fsnow_data, fsnow)
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696 dataloc = 2
01697 sec6hr = secday/c4
01698 maxrec = 1460
01699
01700
01701 recnum = 4*int(yday) - 3 + int(real(sec,kind=dbl_kind)/sec6hr)
01702
01703
01704
01705 ixm = mod(recnum+maxrec-2,maxrec) + 1
01706 ixx = mod(recnum-1, maxrec) + 1
01707
01708
01709
01710
01711
01712
01713 recslot = 2
01714 ixp = 99
01715 call interp_coeff (recnum, recslot, sec6hr, dataloc)
01716
01717
01718 read6 = .false.
01719 if (istep==1 .or. oldrecnum /= recnum) read6 = .true.
01720
01721 if (trim(atm_data_format) == 'bin') then
01722 call read_data (read6, 0, fyear, ixm, ixx, ixp, &
01723 maxrec, tair_file, Tair_data, &
01724 field_loc_center, field_type_scalar)
01725 call read_data (read6, 0, fyear, ixm, ixx, ixp, &
01726 maxrec, uwind_file, uatm_data, &
01727 field_loc_center, field_type_scalar)
01728 call read_data (read6, 0, fyear, ixm, ixx, ixp, &
01729 maxrec, vwind_file, vatm_data, &
01730 field_loc_center, field_type_scalar)
01731 call read_data (read6, 0, fyear, ixm, ixx, ixp, &
01732 maxrec, rhoa_file, rhoa_data, &
01733 field_loc_center, field_type_scalar)
01734 call read_data (read6, 0, fyear, ixm, ixx, ixp, &
01735 maxrec, humid_file, Qa_data, &
01736 field_loc_center, field_type_scalar)
01737 else
01738 call abort_ice ('nonbinary atm_data_format unavailable')
01739 endif
01740
01741
01742 call interpolate_data (Tair_data, Tair)
01743 call interpolate_data (uatm_data, uatm)
01744 call interpolate_data (vatm_data, vatm)
01745 call interpolate_data (rhoa_data, rhoa)
01746 call interpolate_data (Qa_data, Qa)
01747
01748
01749 oldrecnum = recnum
01750
01751 end subroutine ncar_data
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762 subroutine ecmwf_files (yr)
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782 integer (kind=int_kind), intent(in) ::
01783 yr
01784
01785
01786
01787 fsw_file = &
01788 trim(atm_data_dir)//'sol_2002.r'
01789 call file_year(fsw_file,yr)
01790
01791 flw_file = &
01792 trim(atm_data_dir)//'flo_2002.r'
01793 call file_year(flw_file,yr)
01794
01795 rain_file = &
01796 trim(atm_data_dir)//'prec_lanl_12.r'
01797
01798
01799
01800 uwind_file = &
01801 trim(atm_data_dir)//'ucmp_2002.r'
01802 call file_year(uwind_file,yr)
01803
01804 vwind_file = &
01805 trim(atm_data_dir)//'vcmp_2002.r'
01806 call file_year(vwind_file,yr)
01807
01808 tair_file = &
01809 trim(atm_data_dir)//'tair_2002.r'
01810 call file_year(tair_file,yr)
01811
01812 humid_file = &
01813 trim(atm_data_dir)//'qa_2002.r'
01814 call file_year(humid_file,yr)
01815
01816 rhoa_file = &
01817 trim(atm_data_dir)//'rhoa_ncar85-88_12.r'
01818
01819
01820
01821 if (my_task == master_task) then
01822 write (nu_diag,*) ' '
01823 write (nu_diag,*) 'Forcing data year = ', fyear
01824 write (nu_diag,*) 'Atmospheric data files:'
01825 write (nu_diag,*) trim(fsw_file)
01826 write (nu_diag,*) trim(flw_file)
01827 write (nu_diag,*) trim(rain_file)
01828 write (nu_diag,*) trim(uwind_file)
01829 write (nu_diag,*) trim(vwind_file)
01830 write (nu_diag,*) trim(tair_file)
01831 write (nu_diag,*) trim(humid_file)
01832 write (nu_diag,*) trim(rhoa_file)
01833 endif
01834
01835 end subroutine ecmwf_files
01836
01837
01838
01839
01840
01841
01842
01843
01844 subroutine ECMWF_data
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858 use ice_flux
01859
01860
01861
01862
01863
01864
01865 integer (kind=int_kind) ::
01866 ixm,ixx,ixp ,
01867 recnum ,
01868 maxrec ,
01869 recslot ,
01870 dataloc ,
01871
01872 midmonth
01873
01874 logical (kind=log_kind) :: readm, readd
01875
01876
01877
01878
01879
01880
01881
01882
01883 midmonth = 15
01884
01885
01886
01887 maxrec = 12
01888 ixm = mod(month+maxrec-2,maxrec) + 1
01889 ixp = mod(month, maxrec) + 1
01890 if (mday >= midmonth) ixm = 99
01891 if (mday < midmonth) ixp = 99
01892
01893
01894
01895
01896
01897
01898 recslot = 1
01899 if (mday < midmonth) recslot = 2
01900
01901
01902 call interp_coeff_monthly (recslot)
01903
01904
01905 readm = .false.
01906 if (istep==1 .or. (mday==midmonth .and. sec==0)) readm = .true.
01907
01908 if (trim(atm_data_format) == 'bin') then
01909 call read_clim_data (readm, 0, ixm, month, ixp, &
01910 rhoa_file, rhoa_data, field_loc_center, field_type_scalar)
01911 call read_clim_data (readm, 0, ixm, month, ixp, &
01912 rain_file, fsnow_data, field_loc_center, field_type_scalar)
01913 else
01914 call abort_ice ('nonbinary atm_data_format unavailable')
01915 endif
01916
01917
01918 call interpolate_data (fsnow_data, fsnow)
01919 call interpolate_data (rhoa_data, rhoa)
01920
01921
01922
01923
01924
01925
01926
01927
01928 dataloc = 1
01929 maxrec = 365
01930
01931
01932 recnum = int(yday)
01933
01934
01935
01936 ixm = mod(recnum+maxrec-2,maxrec) + 1
01937 ixx = mod(recnum-1, maxrec) + 1
01938 ixp = mod(recnum, maxrec) + 1
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949 if (real(sec,kind=dbl_kind) < p5*secday-puny) then
01950 recslot = 2
01951 ixp = 99
01952 else
01953 recslot = 1
01954 ixm = 99
01955 endif
01956
01957 call interp_coeff (recnum, recslot, secday, dataloc)
01958
01959
01960
01961 readd = .false.
01962 if (istep==1 .or. (recslot==1 .and. oldrecslot==2)) &
01963 readd = .true.
01964
01965 if (trim(atm_data_format) == 'bin') then
01966 call read_data (readd, 0, fyear, ixm, ixx, ixp, maxrec, &
01967 tair_file, Tair_data, &
01968 field_loc_center, field_type_scalar)
01969 call read_data (readd, 0, fyear, ixm, ixx, ixp, maxrec, &
01970 uwind_file, uatm_data, &
01971 field_loc_center, field_type_scalar)
01972 call read_data (readd, 0, fyear, ixm, ixx, ixp, maxrec, &
01973 vwind_file, vatm_data, &
01974 field_loc_center, field_type_scalar)
01975 call read_data (readd, 0, fyear, ixm, ixx, ixp, maxrec, &
01976 fsw_file, fsw_data, &
01977 field_loc_center, field_type_scalar)
01978 call read_data (readd, 0, fyear, ixm, ixx, ixp, maxrec, &
01979 flw_file, flw_data, &
01980 field_loc_center, field_type_scalar)
01981 call read_data (readd, 0, fyear, ixm, ixx, ixp, maxrec, &
01982 humid_file, Qa_data, &
01983 field_loc_center, field_type_scalar)
01984 else
01985 call abort_ice ('nonbinary atm_data_format unavailable')
01986 endif
01987
01988
01989 call interpolate_data (Tair_data, Tair)
01990 call interpolate_data (uatm_data, uatm)
01991 call interpolate_data (vatm_data, vatm)
01992 call interpolate_data ( fsw_data, fsw)
01993 call interpolate_data ( flw_data, flw)
01994 call interpolate_data ( Qa_data, Qa)
01995
01996
01997 oldrecslot = recslot
01998
01999 end subroutine ecmwf_data
02000
02001
02002
02003
02004
02005
02006
02007
02008
02009
02010
02011
02012 subroutine LY_files (yr)
02013
02014
02015
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029
02030
02031 integer (kind=int_kind), intent(in) ::
02032 yr
02033
02034
02035
02036 flw_file = &
02037 trim(atm_data_dir)//'MONTHLY/cldf.omip.dat'
02038
02039 rain_file = &
02040 trim(atm_data_dir)//'MONTHLY/prec.nmyr.dat'
02041
02042 uwind_file = &
02043 trim(atm_data_dir)//'4XDAILY/u_10.1996.dat'
02044 call file_year(uwind_file,yr)
02045
02046 vwind_file = &
02047 trim(atm_data_dir)//'4XDAILY/v_10.1996.dat'
02048 call file_year(vwind_file,yr)
02049
02050 tair_file = &
02051 trim(atm_data_dir)//'4XDAILY/t_10.1996.dat'
02052 call file_year(tair_file,yr)
02053
02054 humid_file = &
02055 trim(atm_data_dir)//'4XDAILY/q_10.1996.dat'
02056 call file_year(humid_file,yr)
02057
02058 if (my_task == master_task) then
02059 write (nu_diag,*) ' '
02060 write (nu_diag,*) 'Forcing data year = ', fyear
02061 write (nu_diag,*) 'Atmospheric data files:'
02062 write (nu_diag,*) trim(flw_file)
02063 write (nu_diag,*) trim(rain_file)
02064 write (nu_diag,*) trim(uwind_file)
02065 write (nu_diag,*) trim(vwind_file)
02066 write (nu_diag,*) trim(tair_file)
02067 write (nu_diag,*) trim(humid_file)
02068 endif
02069
02070 end subroutine LY_files
02071
02072
02073
02074
02075
02076
02077
02078
02079
02080
02081 subroutine LY_data
02082
02083
02084
02085
02086
02087
02088
02089
02090
02091 use ice_global_reductions
02092 use ice_domain, only: nblocks, distrb_info, blocks_ice
02093 use ice_flux
02094 use ice_grid, only: hm, tlon, tlat, tmask, umask
02095
02096
02097
02098
02099
02100 integer (kind=int_kind) ::
02101 i, j ,
02102 imx,ixx,ipx ,
02103 recnum ,
02104 maxrec ,
02105 recslot ,
02106 midmonth ,
02107 dataloc ,
02108
02109 iblk , &
02110 ilo,ihi,jlo,jhi
02111
02112 real (kind=dbl_kind) ::
02113 sec6hr ,
02114 vmin, vmax
02115
02116 logical (kind=log_kind) :: readm, read6
02117
02118 type (block) ::
02119 this_block
02120
02121
02122
02123
02124
02125
02126
02127
02128 midmonth = 15
02129
02130
02131
02132 maxrec = 12
02133 imx = mod(month+maxrec-2,maxrec) + 1
02134 ipx = mod(month, maxrec) + 1
02135 if (mday >= midmonth) imx = 99
02136 if (mday < midmonth) ipx = 99
02137
02138
02139
02140
02141
02142
02143 recslot = 1
02144 if (mday < midmonth) recslot = 2
02145
02146
02147 call interp_coeff_monthly (recslot)
02148
02149
02150 readm = .false.
02151 if (istep==1 .or. (mday==midmonth .and. sec==0)) readm = .true.
02152
02153 call read_clim_data (readm, 0, imx, month, ipx, &
02154 flw_file, cldf_data, field_loc_center, field_type_scalar)
02155 call read_clim_data (readm, 0, imx, month, ipx, &
02156 rain_file, fsnow_data, field_loc_center, field_type_scalar)
02157
02158 call interpolate_data (cldf_data, cldf)
02159 call interpolate_data (fsnow_data, fsnow)
02160
02161
02162
02163
02164
02165
02166
02167
02168
02169 dataloc = 2
02170 sec6hr = secday/c4
02171 maxrec = 1460
02172
02173
02174 recnum = 4*int(yday) - 3 + int(real(sec,kind=dbl_kind)/sec6hr)
02175
02176
02177
02178 imx = mod(recnum+maxrec-2,maxrec) + 1
02179 ixx = mod(recnum-1, maxrec) + 1
02180
02181
02182
02183
02184
02185
02186 recslot = 2
02187 ipx = 99
02188 call interp_coeff (recnum, recslot, sec6hr, dataloc)
02189
02190
02191 read6 = .false.
02192 if (istep==1 .or. oldrecnum .ne. recnum) read6 = .true.
02193
02194 if (trim(atm_data_format) == 'bin') then
02195 call read_data (read6, 0, fyear, imx, ixx, ipx, maxrec, &
02196 tair_file, Tair_data, &
02197 field_loc_center, field_type_scalar)
02198 call read_data (read6, 0, fyear, imx, ixx, ipx, maxrec, &
02199 uwind_file, uatm_data, &
02200 field_loc_center, field_type_scalar)
02201 call read_data (read6, 0, fyear, imx, ixx, ipx, maxrec, &
02202 vwind_file, vatm_data, &
02203 field_loc_center, field_type_scalar)
02204 call read_data (read6, 0, fyear, imx, ixx, ipx, maxrec, &
02205 humid_file, Qa_data, &
02206 field_loc_center, field_type_scalar)
02207 else
02208 call abort_ice ('nonbinary atm_data_format unavailable')
02209 endif
02210
02211
02212 call interpolate_data (Tair_data, Tair)
02213 call interpolate_data (uatm_data, uatm)
02214 call interpolate_data (vatm_data, vatm)
02215 call interpolate_data (Qa_data, Qa)
02216
02217 do iblk = 1, nblocks
02218 call Qa_fixLY(nx_block, ny_block, &
02219 Tair (:,:,iblk), &
02220 Qa (:,:,iblk))
02221
02222 do j = 1, ny_block
02223 do i = 1, nx_block
02224 Qa (i,j,iblk) = Qa (i,j,iblk) * hm(i,j,iblk)
02225 Tair(i,j,iblk) = Tair(i,j,iblk) * hm(i,j,iblk)
02226 uatm(i,j,iblk) = uatm(i,j,iblk) * hm(i,j,iblk)
02227 vatm(i,j,iblk) = vatm(i,j,iblk) * hm(i,j,iblk)
02228 enddo
02229 enddo
02230
02231
02232 this_block = get_block(blocks_ice(iblk),iblk)
02233 ilo = this_block%ilo
02234 ihi = this_block%ihi
02235 jlo = this_block%jlo
02236 jhi = this_block%jhi
02237
02238 call compute_shortwave(nx_block, ny_block, &
02239 ilo, ihi, jlo, jhi, &
02240 TLON (:,:,iblk), &
02241 TLAT (:,:,iblk), &
02242 hm (:,:,iblk), &
02243 Qa (:,:,iblk), &
02244 cldf (:,:,iblk), &
02245 fsw (:,:,iblk))
02246
02247 enddo
02248
02249
02250 oldrecnum = recnum
02251
02252 if (dbug) then
02253 if (my_task == master_task) write (nu_diag,*) 'LY_bulk_data'
02254 vmin = global_minval(fsw,distrb_info,tmask)
02255
02256 vmax = global_maxval(fsw,distrb_info,tmask)
02257 if (my_task.eq.master_task) &
02258 write (nu_diag,*) 'fsw',vmin,vmax
02259 vmin = global_minval(cldf,distrb_info,tmask)
02260 vmax = global_maxval(cldf,distrb_info,tmask)
02261 if (my_task.eq.master_task) &
02262 write (nu_diag,*) 'cldf',vmin,vmax
02263 vmin =global_minval(fsnow,distrb_info,tmask)
02264 vmax =global_maxval(fsnow,distrb_info,tmask)
02265 if (my_task.eq.master_task) &
02266 write (nu_diag,*) 'fsnow',vmin,vmax
02267 vmin = global_minval(Tair,distrb_info,tmask)
02268 vmax = global_maxval(Tair,distrb_info,tmask)
02269 if (my_task.eq.master_task) &
02270 write (nu_diag,*) 'Tair',vmin,vmax
02271 vmin = global_minval(uatm,distrb_info,umask)
02272 vmax = global_maxval(uatm,distrb_info,umask)
02273 if (my_task.eq.master_task) &
02274 write (nu_diag,*) 'uatm',vmin,vmax
02275 vmin = global_minval(vatm,distrb_info,umask)
02276 vmax = global_maxval(vatm,distrb_info,umask)
02277 if (my_task.eq.master_task) &
02278 write (nu_diag,*) 'vatm',vmin,vmax
02279 vmin = global_minval(Qa,distrb_info,tmask)
02280 vmax = global_maxval(Qa,distrb_info,tmask)
02281 if (my_task.eq.master_task) &
02282 write (nu_diag,*) 'Qa',vmin,vmax
02283
02284 endif
02285
02286 end subroutine LY_data
02287
02288
02289
02290 subroutine compute_shortwave(nx_block, ny_block, &
02291 ilo, ihi, jlo, jhi, &
02292 TLON, TLAT, hm, Qa, cldf, fsw)
02293
02294
02295
02296
02297
02298 integer (kind=int_kind), intent(in) ::
02299 nx_block, ny_block,
02300 ilo,ihi,jlo,jhi
02301
02302 real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) ::
02303 TLON, TLAT ,
02304 Qa ,
02305 cldf ,
02306 hm
02307
02308 real (kind=dbl_kind), dimension(nx_block,ny_block),
02309 intent(inout) ::
02310 fsw
02311
02312 real (kind=dbl_kind) ::
02313 hour_angle,
02314 solar_time,
02315 declin ,
02316 cosZ ,
02317 e, d ,
02318 sw0 ,
02319 deg2rad
02320
02321 integer (kind=int_kind) ::
02322 i, j
02323
02324 do j=jlo,jhi
02325 do i=ilo,ihi
02326 deg2rad = pi/c180
02327 solar_time = mod(real(sec,kind=dbl_kind),secday)/c3600
02328 + c12*sin(p5*TLON(i,j))
02329 hour_angle = (c12 - solar_time)*pi/c12
02330 declin = 23.44_dbl_kind*cos((172._dbl_kind-yday) &
02331 * c2*pi/c365)*deg2rad
02332 cosZ = sin(TLAT(i,j))*sin(declin) &
02333 + cos(TLAT(i,j))*cos(declin)*cos(hour_angle)
02334 cosZ = max(cosZ,c0)
02335 e = 1.e5*Qa(i,j)/(0.622_dbl_kind + 0.378_dbl_kind*Qa(i,j))
02336 d = (cosZ+2.7_dbl_kind)*e*1.e-5_dbl_kind+1.085_dbl_kind*cosZ+p1
02337 sw0 = 1353._dbl_kind*cosZ**2/d
02338 sw0 = max(sw0,c0)
02339
02340
02341 Fsw(i,j) = sw0*(c1-p6*cldf(i,j)**3)
02342 Fsw(i,j) = Fsw(i,j)*hm(i,j)
02343 enddo
02344 enddo
02345
02346 end subroutine compute_shortwave
02347
02348
02349
02350 subroutine Qa_fixLY(nx_block, ny_block, Tair, Qa)
02351
02352 use ice_work, only: worka
02353
02354 integer (kind=int_kind), intent(in) ::
02355 nx_block, ny_block
02356
02357 real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) ::
02358 Tair
02359
02360 real (kind=dbl_kind), dimension(nx_block,ny_block),
02361 intent(inout) ::
02362 Qa
02363
02364 worka = Tair - Tffresh
02365 worka = c2 + (0.7859_dbl_kind + 0.03477_dbl_kind*worka) &
02366 /(c1 + 0.00412_dbl_kind*worka) &
02367 + 0.00422_dbl_kind*worka
02368
02369 worka = (c10**worka)
02370 worka = max(worka,puny)
02371
02372 worka = 0.622_dbl_kind*worka/(1.e5_dbl_kind-0.378_dbl_kind*worka)
02373
02374 Qa = min(Qa, worka)
02375
02376 end subroutine Qa_fixLY
02377
02378
02379
02380
02381
02382
02383
02384
02385
02386
02387
02388 subroutine hadgem_files (yr)
02389
02390
02391
02392
02393
02394
02395
02396
02397
02398
02399
02400
02401
02402
02403 use ice_therm_vertical, only: calc_Tsfc
02404 use ice_ocean, only: oceanmixed_ice
02405
02406
02407
02408 integer (kind=int_kind), intent(in) ::
02409 yr
02410
02411
02412
02413 integer (kind=int_kind) ::
02414 n
02415
02416
02417
02418
02419
02420 snow_file = &
02421 trim(atm_data_dir)//'MONTHLY/snowfall.1996.nc'
02422 call file_year(snow_file,yr)
02423
02424 rain_file = &
02425 trim(atm_data_dir)//'MONTHLY/rainfall.1996.nc'
02426 call file_year(rain_file,yr)
02427
02428 if (my_task == master_task) then
02429 write (nu_diag,*) ' '
02430 write (nu_diag,*) 'Atmospheric data files:'
02431 write (nu_diag,*) trim(rain_file)
02432 write (nu_diag,*) trim(snow_file)
02433 endif
02434
02435 if (calc_strair) then
02436
02437
02438
02439
02440
02441 uwind_file = &
02442 trim(atm_data_dir)//'MONTHLY/u_10.1996.nc'
02443 call file_year(uwind_file,yr)
02444
02445 vwind_file = &
02446 trim(atm_data_dir)//'MONTHLY/v_10.1996.nc'
02447 call file_year(vwind_file,yr)
02448
02449 if (my_task == master_task) then
02450 write (nu_diag,*) trim(uwind_file)
02451 write (nu_diag,*) trim(vwind_file)
02452 endif
02453
02454 else
02455
02456
02457
02458
02459
02460 strax_file = &
02461 trim(atm_data_dir)//'MONTHLY/taux.1996.nc'
02462 call file_year(strax_file,yr)
02463
02464 stray_file = &
02465 trim(atm_data_dir)//'MONTHLY/tauy.1996.nc'
02466 call file_year(stray_file,yr)
02467
02468 if (my_task == master_task) then
02469 write (nu_diag,*) trim(strax_file)
02470 write (nu_diag,*) trim(stray_file)
02471 endif
02472
02473 if (calc_Tsfc .or. oceanmixed_ice) then
02474
02475
02476
02477
02478
02479 wind_file = &
02480 trim(atm_data_dir)//'MONTHLY/wind_10.1996.nc'
02481 call file_year(wind_file,yr)
02482
02483 if (my_task == master_task) then
02484 write (nu_diag,*) trim(wind_file)
02485 endif
02486
02487 endif
02488
02489 endif
02490
02491
02492
02493
02494
02495
02496
02497
02498 if (calc_Tsfc .or. oceanmixed_ice .or. calc_strair) then
02499
02500 fsw_file = &
02501 trim(atm_data_dir)//'MONTHLY/SW_incoming.1996.nc'
02502 call file_year(fsw_file,yr)
02503
02504 flw_file = &
02505 trim(atm_data_dir)//'MONTHLY/LW_incoming.1996.nc'
02506 call file_year(flw_file,yr)
02507
02508 tair_file = &
02509 trim(atm_data_dir)//'MONTHLY/t_10.1996.nc'
02510 call file_year(tair_file,yr)
02511
02512 humid_file = &
02513 trim(atm_data_dir)//'MONTHLY/q_10.1996.nc'
02514 call file_year(humid_file,yr)
02515
02516 rhoa_file = &
02517 trim(atm_data_dir)//'MONTHLY/rho_10.1996.nc'
02518 call file_year(rhoa_file,yr)
02519
02520 if (my_task == master_task) then
02521 write (nu_diag,*) trim(fsw_file)
02522 write (nu_diag,*) trim(flw_file)
02523 write (nu_diag,*) trim(tair_file)
02524 write (nu_diag,*) trim(humid_file)
02525 write (nu_diag,*) trim(rhoa_file)
02526 endif
02527
02528 endif
02529
02530 if (.not. calc_Tsfc) then
02531
02532
02533
02534
02535
02536 do n = 1, ncat
02537
02538
02539 write(topmelt_file(n), '(a,i1,a)') &
02540 trim(atm_data_dir)//'MONTHLY/topmeltn',n,'.1996.nc'
02541 call file_year(topmelt_file(n),yr)
02542
02543
02544 write(botmelt_file(n), '(a,i1,a)') &
02545 trim(atm_data_dir)//'MONTHLY/botmeltn',n,'.1996.nc'
02546 call file_year(botmelt_file(n),yr)
02547
02548 enddo
02549
02550
02551 sublim_file = &
02552 trim(atm_data_dir)//'MONTHLY/sublim.1996.nc'
02553 call file_year(sublim_file,yr)
02554
02555 if (my_task == master_task) then
02556 do n = 1, ncat
02557 write (nu_diag,*) trim(topmelt_file(n))
02558 write (nu_diag,*) trim(botmelt_file(n))
02559 enddo
02560 write (nu_diag,*) trim(sublim_file)
02561
02562 endif
02563
02564 endif
02565
02566 end subroutine hadgem_files
02567
02568
02569
02570
02571
02572
02573
02574
02575
02576 subroutine hadgem_data
02577
02578
02579
02580
02581
02582
02583
02584
02585
02586 use ice_domain, only: nblocks
02587 use ice_flux
02588 use ice_state, only: aice,aicen
02589 use ice_ocean, only: oceanmixed_ice
02590 use ice_therm_vertical, only: calc_Tsfc
02591
02592
02593
02594
02595
02596
02597 integer (kind=int_kind) ::
02598 i, j ,
02599 n ,
02600 iblk ,
02601 ixm,ixx,ixp ,
02602 recnum ,
02603 maxrec ,
02604 recslot ,
02605 dataloc ,
02606
02607 midmonth
02608
02609 logical (kind=log_kind) :: readm
02610
02611 real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) ::
02612 topmelt,
02613 botmelt,
02614 sublim
02615
02616 character (char_len) ::
02617 fieldname
02618
02619
02620
02621
02622
02623
02624
02625
02626 midmonth = 15
02627
02628
02629
02630 maxrec = 12
02631 ixm = mod(month+maxrec-2,maxrec) + 1
02632 ixp = mod(month, maxrec) + 1
02633 if (mday >= midmonth) ixm = 99
02634 if (mday < midmonth) ixp = 99
02635
02636
02637
02638
02639
02640
02641 recslot = 1
02642 if (mday < midmonth) recslot = 2
02643
02644
02645 call interp_coeff_monthly (recslot)
02646
02647
02648 readm = .false.
02649 if (istep==1 .or. (mday==midmonth .and. sec==0)) readm = .true.
02650
02651
02652
02653
02654
02655 fieldname='rainfall'
02656 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
02657 maxrec, rain_file, fieldname, frain_data, &
02658 field_loc_center, field_type_scalar)
02659 fieldname='snowfall'
02660 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
02661 maxrec, snow_file, fieldname, fsnow_data, &
02662 field_loc_center, field_type_scalar)
02663
02664
02665 call interpolate_data (fsnow_data, fsnow)
02666 call interpolate_data (frain_data, frain)
02667
02668 if (calc_strair) then
02669
02670
02671
02672
02673
02674 fieldname='u_10'
02675 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
02676 maxrec, uwind_file, fieldname, uatm_data, &
02677 field_loc_center, field_type_vector)
02678 fieldname='v_10'
02679 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
02680 maxrec, vwind_file, fieldname, vatm_data, &
02681 field_loc_center, field_type_vector)
02682
02683
02684 call interpolate_data (uatm_data, uatm)
02685 call interpolate_data (vatm_data, vatm)
02686
02687 else
02688
02689
02690
02691
02692
02693 fieldname='taux'
02694 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
02695 maxrec, strax_file, fieldname, strax_data, &
02696 field_loc_center, field_type_vector)
02697 fieldname='tauy'
02698 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
02699 maxrec, stray_file, fieldname, stray_data, &
02700 field_loc_center, field_type_vector)
02701
02702
02703 call interpolate_data (strax_data, strax)
02704 call interpolate_data (stray_data, stray)
02705
02706 if (calc_Tsfc .or. oceanmixed_ice) then
02707
02708
02709
02710
02711
02712 fieldname='wind_10'
02713 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
02714 maxrec, wind_file, fieldname, wind_data, &
02715 field_loc_center, field_type_scalar)
02716
02717
02718 call interpolate_data (wind_data, wind)
02719
02720 endif
02721
02722 endif
02723
02724
02725
02726
02727
02728
02729
02730
02731
02732
02733 if (calc_Tsfc .or. oceanmixed_ice .or. calc_strair) then
02734
02735 fieldname='SW_incoming'
02736 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
02737 maxrec, fsw_file, fieldname, fsw_data, &
02738 field_loc_center, field_type_scalar)
02739 fieldname='LW_incoming'
02740 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
02741 maxrec, flw_file, fieldname, flw_data, &
02742 field_loc_center, field_type_scalar)
02743 fieldname='t_10'
02744 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
02745 maxrec, tair_file, fieldname, Tair_data, &
02746 field_loc_center, field_type_scalar)
02747 fieldname='rho_10'
02748 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
02749 maxrec, rhoa_file, fieldname, rhoa_data, &
02750 field_loc_center, field_type_scalar)
02751 fieldname='q_10'
02752 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
02753 maxrec, humid_file, fieldname, Qa_data, &
02754 field_loc_center, field_type_scalar)
02755
02756
02757
02758 call interpolate_data (fsw_data, fsw)
02759 call interpolate_data (flw_data, flw)
02760 call interpolate_data (Tair_data, Tair)
02761 call interpolate_data (rhoa_data, rhoa)
02762 call interpolate_data (Qa_data, Qa)
02763
02764 endif
02765
02766 if (.not. calc_Tsfc) then
02767
02768
02769
02770
02771
02772 fieldname='sublim'
02773 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
02774 maxrec, sublim_file, fieldname, sublim_data, &
02775 field_loc_center, field_type_scalar)
02776
02777
02778 call interpolate_data (sublim_data, sublim)
02779
02780 do n = 1, ncat
02781 write(fieldname, '(a,i1)') 'topmeltn',n
02782 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
02783 maxrec, topmelt_file(n), fieldname, topmelt_data(:,:,:,:,n), &
02784 field_loc_center, field_type_scalar)
02785
02786 write(fieldname, '(a,i1)') 'botmeltn',n
02787 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
02788 maxrec, botmelt_file(n), fieldname, botmelt_data(:,:,:,:,n), &
02789 field_loc_center, field_type_scalar)
02790
02791 call interpolate_data (topmelt_data(:,:,:,:,n), topmelt)
02792 call interpolate_data (botmelt_data(:,:,:,:,n), botmelt)
02793
02794
02795
02796
02797
02798
02799
02800
02801
02802
02803
02804 do iblk = 1, nblocks
02805 do j = 1, ny_block
02806 do i = 1, nx_block
02807
02808 fcondtopn_f(i,j,n,iblk) = botmelt(i,j,iblk)
02809 fsurfn_f(i,j,n,iblk) = topmelt(i,j,iblk) &
02810 + botmelt(i,j,iblk)
02811 flatn_f(i,j,n,iblk) = - sublim(i,j,iblk)*Lsub
02812
02813 enddo
02814 enddo
02815
02816 enddo
02817
02818 enddo
02819
02820 endif
02821
02822 end subroutine hadgem_data
02823
02824
02825
02826
02827
02828
02829
02830
02831
02832
02833
02834 subroutine monthly_files (yr)
02835
02836
02837
02838
02839
02840
02841
02842
02843
02844
02845
02846
02847
02848
02849
02850
02851
02852 integer (kind=int_kind), intent(in) ::
02853 yr
02854
02855
02856
02857 flw_file = &
02858 trim(atm_data_dir)//'MONTHLY/cldf.omip.dat'
02859
02860 rain_file = &
02861 trim(atm_data_dir)//'MONTHLY/prec.nmyr.dat'
02862
02863 tair_file = &
02864 trim(atm_data_dir)//'MONTHLY/t_10.1996.dat'
02865 call file_year(tair_file,yr)
02866
02867 humid_file = &
02868 trim(atm_data_dir)//'MONTHLY/q_10.1996.dat'
02869 call file_year(humid_file,yr)
02870
02871
02872 strax_file = &
02873 trim(atm_data_dir)//'MONTHLY/strx.1996.dat'
02874 call file_year(strax_file,yr)
02875
02876 stray_file = &
02877 trim(atm_data_dir)//'MONTHLY/stry.1996.dat'
02878 call file_year(stray_file,yr)
02879
02880 wind_file = &
02881 trim(atm_data_dir)//'MONTHLY/wind.1996.dat'
02882 call file_year(wind_file,yr)
02883
02884 if (my_task == master_task) then
02885 write (nu_diag,*) ' '
02886 write (nu_diag,*) 'Forcing data year = ', fyear
02887 write (nu_diag,*) 'Atmospheric data files:'
02888 write (nu_diag,*) trim(flw_file)
02889 write (nu_diag,*) trim(rain_file)
02890 write (nu_diag,*) trim(tair_file)
02891 write (nu_diag,*) trim(humid_file)
02892 write (nu_diag,*) trim(uwind_file)
02893 write (nu_diag,*) trim(vwind_file)
02894 endif
02895
02896 end subroutine monthly_files
02897
02898
02899
02900
02901
02902
02903
02904
02905
02906 subroutine monthly_data
02907
02908
02909
02910
02911
02912
02913
02914
02915
02916 use ice_global_reductions
02917 use ice_domain, only: nblocks, distrb_info, blocks_ice
02918 use ice_flux
02919 use ice_grid, only: hm, tlon, tlat, tmask, umask
02920
02921
02922
02923
02924
02925 integer (kind=int_kind) ::
02926 i, j ,
02927 imx,ipx ,
02928 recnum ,
02929 maxrec ,
02930 recslot ,
02931 midmonth ,
02932 iblk ,
02933 ilo,ihi,jlo,jhi
02934
02935 real (kind=dbl_kind) ::
02936 vmin, vmax
02937
02938 logical (kind=log_kind) :: readm
02939
02940 type (block) ::
02941 this_block
02942
02943
02944
02945
02946
02947
02948
02949
02950 midmonth = 15
02951
02952
02953
02954 maxrec = 12
02955 imx = mod(month+maxrec-2,maxrec) + 1
02956 ipx = mod(month, maxrec) + 1
02957 if (mday >= midmonth) imx = 99
02958 if (mday < midmonth) ipx = 99
02959
02960
02961
02962
02963
02964
02965 recslot = 1
02966 if (mday < midmonth) recslot = 2
02967
02968
02969 call interp_coeff_monthly (recslot)
02970
02971
02972 readm = .false.
02973 if (istep==1 .or. (mday==midmonth .and. sec==0)) readm = .true.
02974
02975 call read_clim_data (readm, 0, imx, month, ipx, &
02976 flw_file, cldf_data, &
02977 field_loc_center, field_type_scalar)
02978 call read_clim_data (readm, 0, imx, month, ipx, &
02979 rain_file, fsnow_data, &
02980 field_loc_center, field_type_scalar)
02981 call read_clim_data (readm, 0, imx, month, ipx, &
02982 tair_file, Tair_data, &
02983 field_loc_center, field_type_scalar)
02984 call read_clim_data (readm, 0, imx, month, ipx, &
02985 humid_file, Qa_data, &
02986 field_loc_center, field_type_scalar)
02987 call read_clim_data (readm, 0, imx, month, ipx, &
02988 wind_file, wind_data, &
02989 field_loc_center, field_type_scalar)
02990 call read_clim_data (readm, 0, imx, month, ipx, &
02991 strax_file, strax_data, &
02992 field_loc_center, field_type_vector)
02993 call read_clim_data (readm, 0, imx, month, ipx, &
02994 stray_file, stray_data, &
02995 field_loc_center, field_type_vector)
02996
02997 call interpolate_data (cldf_data, cldf)
02998 call interpolate_data (fsnow_data, fsnow)
02999 call interpolate_data (Tair_data, Tair)
03000 call interpolate_data (Qa_data, Qa)
03001 call interpolate_data (wind_data, wind)
03002 call interpolate_data (strax_data, strax)
03003 call interpolate_data (stray_data, stray)
03004
03005 do iblk = 1, nblocks
03006 call Qa_fixLY(nx_block, ny_block, &
03007 Tair (:,:,iblk), &
03008 Qa (:,:,iblk))
03009
03010 do j = 1, ny_block
03011 do i = 1, nx_block
03012 Qa (i,j,iblk) = Qa (i,j,iblk) * hm(i,j,iblk)
03013 Tair (i,j,iblk) = Tair (i,j,iblk) * hm(i,j,iblk)
03014 wind (i,j,iblk) = wind (i,j,iblk) * hm(i,j,iblk)
03015 strax(i,j,iblk) = strax(i,j,iblk) * hm(i,j,iblk)
03016 stray(i,j,iblk) = stray(i,j,iblk) * hm(i,j,iblk)
03017 enddo
03018 enddo
03019
03020
03021 this_block = get_block(blocks_ice(iblk),iblk)
03022 ilo = this_block%ilo
03023 ihi = this_block%ihi
03024 jlo = this_block%jlo
03025 jhi = this_block%jhi
03026
03027 call compute_shortwave(nx_block, ny_block, &
03028 ilo, ihi, jlo, jhi, &
03029 TLON (:,:,iblk), &
03030 TLAT (:,:,iblk), &
03031 hm (:,:,iblk), &
03032 Qa (:,:,iblk), &
03033 cldf (:,:,iblk), &
03034 fsw (:,:,iblk))
03035
03036 enddo
03037
03038
03039 oldrecnum = recnum
03040
03041 if (dbug) then
03042 if (my_task == master_task) write (nu_diag,*) 'LY_bulk_data'
03043 vmin = global_minval(fsw,distrb_info,tmask)
03044 vmax = global_maxval(fsw,distrb_info,tmask)
03045 if (my_task.eq.master_task) &
03046 write (nu_diag,*) 'fsw',vmin,vmax
03047 vmin = global_minval(cldf,distrb_info,tmask)
03048 vmax = global_maxval(cldf,distrb_info,tmask)
03049 if (my_task.eq.master_task) &
03050 write (nu_diag,*) 'cldf',vmin,vmax
03051 vmin =global_minval(fsnow,distrb_info,tmask)
03052 vmax =global_maxval(fsnow,distrb_info,tmask)
03053 if (my_task.eq.master_task) &
03054 write (nu_diag,*) 'fsnow',vmin,vmax
03055 vmin = global_minval(Tair,distrb_info,tmask)
03056 vmax = global_maxval(Tair,distrb_info,tmask)
03057 if (my_task.eq.master_task) &
03058 write (nu_diag,*) 'Tair',vmin,vmax
03059 vmin = global_minval(wind,distrb_info,umask)
03060 vmax = global_maxval(wind,distrb_info,umask)
03061 if (my_task.eq.master_task) &
03062 write (nu_diag,*) 'wind',vmin,vmax
03063 vmin = global_minval(strax,distrb_info,umask)
03064 vmax = global_maxval(strax,distrb_info,umask)
03065 if (my_task.eq.master_task) &
03066 write (nu_diag,*) 'strax',vmin,vmax
03067 vmin = global_minval(stray,distrb_info,umask)
03068 vmax = global_maxval(stray,distrb_info,umask)
03069 if (my_task.eq.master_task) &
03070 write (nu_diag,*) 'stray',vmin,vmax
03071 vmin = global_minval(Qa,distrb_info,tmask)
03072 vmax = global_maxval(Qa,distrb_info,tmask)
03073 if (my_task.eq.master_task) &
03074 write (nu_diag,*) 'Qa',vmin,vmax
03075
03076 endif
03077
03078 end subroutine monthly_data
03079
03080
03081
03082
03083
03084
03085
03086
03087
03088
03089 subroutine ocn_data_clim (dt)
03090
03091
03092
03093
03094
03095
03096
03097
03098
03099
03100
03101
03102
03103 use ice_domain, only: nblocks
03104 use ice_ocean
03105 use ice_flux, only: Tf, sss, sst, uocn, vocn, ss_tltx, ss_tlty, Tfrzpt
03106 use ice_grid, only: t2ugrid_vector
03107
03108
03109
03110 real (kind=dbl_kind), intent(in) ::
03111 dt
03112
03113
03114
03115 integer (kind=int_kind) ::
03116 i, j, iblk ,
03117 ixm,ixp ,
03118 maxrec ,
03119 recslot ,
03120 midmonth
03121
03122 real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) ::
03123 sstdat
03124
03125 logical (kind=log_kind) :: readm
03126
03127 if (my_task == master_task .and. istep == 1) then
03128 if (trim(sss_data_type)=='clim') then
03129 write (nu_diag,*) ' '
03130 write (nu_diag,*) 'SSS data interpolated to timestep:'
03131 write (nu_diag,*) trim(sss_file)
03132 endif
03133 if (trim(sst_data_type)=='clim') then
03134 write (nu_diag,*) ' '
03135 write (nu_diag,*) 'SST data interpolated to timestep:'
03136 write (nu_diag,*) trim(sst_file)
03137 if (restore_sst) write (nu_diag,*) &
03138 'SST restoring timescale (days) =', trestore
03139 endif
03140 endif
03141
03142
03143
03144
03145
03146
03147
03148
03149 if (trim(sss_data_type)=='clim' .or. &
03150 trim(sst_data_type)=='clim') then
03151
03152 midmonth = 15
03153
03154
03155
03156 maxrec = 12
03157 ixm = mod(month+maxrec-2,maxrec) + 1
03158 ixp = mod(month, maxrec) + 1
03159 if (mday >= midmonth) ixm = 99
03160 if (mday < midmonth) ixp = 99
03161
03162
03163
03164
03165
03166
03167 recslot = 1
03168 if (mday < midmonth) recslot = 2
03169
03170
03171 call interp_coeff_monthly (recslot)
03172
03173 readm = .false.
03174 if (istep==1 .or. (mday==midmonth .and. sec==0)) readm = .true.
03175
03176 endif
03177
03178
03179
03180
03181
03182
03183 if (trim(sss_data_type)=='clim') then
03184 call read_clim_data (readm, 0, ixm, month, ixp, &
03185 sss_file, sss_data, &
03186 field_loc_center, field_type_scalar)
03187 call interpolate_data (sss_data, sss)
03188
03189 do iblk = 1, nblocks
03190 do j = 1, ny_block
03191 do i = 1, nx_block
03192 sss(i,j,iblk) = max(sss(i,j,iblk), c0)
03193 if (trim(Tfrzpt) == 'constant') then
03194 Tf (i,j,iblk) = -1.8_dbl_kind
03195 else
03196 Tf (i,j,iblk) = -depressT * sss(i,j,iblk)
03197 endif
03198 enddo
03199 enddo
03200 enddo
03201 endif
03202
03203
03204
03205
03206
03207
03208 if (trim(sst_data_type)=='clim') then
03209 call read_clim_data (readm, 0, ixm, month, ixp, &
03210 sst_file, sst_data, &
03211 field_loc_center, field_type_scalar)
03212 call interpolate_data (sst_data, sstdat)
03213
03214 if (restore_sst) then
03215 do iblk = 1, nblocks
03216 do j = 1, ny_block
03217 do i = 1, nx_block
03218 sst(i,j,iblk) = sst(i,j,iblk) &
03219 + (sstdat(i,j,iblk)-sst(i,j,iblk))*dt/trest
03220 enddo
03221 enddo
03222 enddo
03223 endif
03224 endif
03225
03226 end subroutine ocn_data_clim
03227
03228
03229
03230
03231
03232
03233
03234
03235
03236
03237
03238 subroutine ocn_data_ncar_init
03239
03240
03241
03242
03243
03244
03245
03246
03247
03248
03249
03250
03251
03252
03253
03254
03255
03256
03257
03258
03259
03260
03261
03262
03263
03264
03265
03266
03267
03268 use ice_domain, only: nblocks, distrb_info
03269 use ice_gather_scatter
03270 use ice_exit
03271 use ice_work, only: work1
03272 use ice_read_write
03273 #ifdef ncdf
03274 use netcdf
03275 #endif
03276
03277
03278
03279
03280
03281 integer (kind=int_kind) ::
03282 n ,
03283 m ,
03284 nrec,
03285 nbits
03286
03287 character(char_len) ::
03288 vname(nfld)
03289 data vname / &
03290 'T', 'S', 'hblt', 'U', 'V', &
03291 'dhdx', 'dhdy', 'qdp' /
03292
03293 integer (kind=int_kind) ::
03294 fid ,
03295 dimid
03296
03297 integer (kind=int_kind) ::
03298 status ,
03299 nlat ,
03300 nlon
03301
03302 if (my_task == master_task) then
03303
03304 write (nu_diag,*) 'WARNING: evp_prep calculates surface tilt'
03305 write (nu_diag,*) 'WARNING: stress from geostrophic currents,'
03306 write (nu_diag,*) 'WARNING: not data from ocean forcing file.'
03307 write (nu_diag,*) 'WARNING: Alter ice_dyn_evp.F if desired.'
03308
03309 if (restore_sst) write (nu_diag,*) &
03310 'SST restoring timescale = ',trestore,' days'
03311
03312 sst_file = trim(ocn_data_dir)//oceanmixed_file
03313
03314
03315
03316
03317 write (nu_diag,*) 'ocean mixed layer forcing data file = ', &
03318 sst_file
03319
03320 endif
03321
03322 if (trim(ocn_data_format) == 'nc') then
03323 #ifdef ncdf
03324 if (my_task == master_task) then
03325 call ice_open_nc(sst_file, fid)
03326
03327
03328 status = nf90_inq_dimid(fid,'ni',dimid)
03329 status = nf90_inquire_dimension(fid,dimid,len=nlon)
03330
03331
03332 status = nf90_inq_dimid(fid,'nj',dimid)
03333 status = nf90_inquire_dimension(fid,dimid,len=nlat)
03334
03335 if( nlon .ne. nx_global ) then
03336 call abort_ice ('ice: ocn frc file nlon ne nx_global')
03337 endif
03338 if( nlat .ne. ny_global ) then
03339 call abort_ice ('ice: ocn frc file nlat ne ny_global')
03340 endif
03341
03342 endif
03343
03344
03345 do n=1,nfld
03346 do m=1,12
03347
03348
03349 if (n >= 4 .and. n <= 7) then
03350 call ice_read_nc(fid, m, vname(n), work1, dbug, &
03351 field_loc_NEcorner, field_type_vector)
03352 else
03353 call ice_read_nc(fid, m, vname(n), work1, dbug, &
03354 field_loc_center, field_type_scalar)
03355 endif
03356 ocn_frc_m(:,:,:,n,m) = work1(:,:,:)
03357
03358 enddo
03359 enddo
03360
03361 if (my_task == master_task) status = nf90_close(fid)
03362 #endif
03363
03364 else
03365
03366 nbits = 64
03367 call ice_open (nu_forcing, sst_file, nbits)
03368
03369 nrec = 0
03370 do n=1,nfld
03371 do m=1,12
03372 nrec = nrec + 1
03373 if (n >= 4 .and. n <= 7) then
03374 call ice_read (nu_forcing, nrec, work1, 'rda8', dbug, &
03375 field_loc_NEcorner, field_type_vector)
03376 else
03377 call ice_read (nu_forcing, nrec, work1, 'rda8', dbug, &
03378 field_loc_center, field_type_scalar)
03379 endif
03380 ocn_frc_m(:,:,:,n,m) = work1(:,:,:)
03381 enddo
03382 enddo
03383 close (nu_forcing)
03384
03385 endif
03386
03387
03388 ocn_frc_m(:,:,:,4,:) = c0
03389 ocn_frc_m(:,:,:,5,:) = c0
03390
03391
03392 end subroutine ocn_data_ncar_init
03393
03394
03395
03396
03397
03398
03399
03400
03401
03402 subroutine ocn_data_ncar(dt)
03403
03404
03405
03406
03407
03408
03409
03410
03411
03412
03413
03414
03415 use ice_global_reductions
03416 use ice_domain, only: nblocks, distrb_info
03417 use ice_gather_scatter
03418 use ice_exit
03419 use ice_work, only: work_g1, work1
03420 use ice_flux, only: sss, sst, Tf, uocn, vocn, ss_tltx, ss_tlty, &
03421 qdp, hmix, Tfrzpt
03422
03423 use ice_restart, only: restart
03424 use ice_grid, only: hm, tmask, umask
03425
03426
03427
03428
03429 real (kind=dbl_kind), intent(in) ::
03430 dt
03431
03432
03433
03434 integer (kind=int_kind) ::
03435 i, j, n, iblk ,
03436 imx,ipx ,
03437 maxrec ,
03438 recslot ,
03439 midmonth
03440
03441 real (kind=dbl_kind) ::
03442 vmin, vmax
03443
03444
03445
03446
03447
03448
03449
03450
03451 midmonth = 15
03452
03453
03454
03455 maxrec = 12
03456 imx = mod(month+maxrec-2,maxrec) + 1
03457 ipx = mod(month, maxrec) + 1
03458 if (mday >= midmonth) imx = 99
03459 if (mday < midmonth) ipx = 99
03460
03461
03462
03463
03464
03465
03466 recslot = 1
03467 if (mday < midmonth) recslot = 2
03468
03469
03470 call interp_coeff_monthly (recslot)
03471
03472 do n = nfld, 1, -1
03473 do iblk = 1, nblocks
03474
03475 if (imx /= 99) then
03476 sst_data(:,:,1,iblk) = ocn_frc_m(:,:,iblk,n,imx)
03477 sst_data(:,:,2,iblk) = ocn_frc_m(:,:,iblk,n,month)
03478 else
03479 sst_data(:,:,1,iblk) = ocn_frc_m(:,:,iblk,n,month)
03480 sst_data(:,:,2,iblk) = ocn_frc_m(:,:,iblk,n,ipx)
03481 endif
03482 enddo
03483 call interpolate_data (sst_data,work1)
03484
03485 do j = 1, ny_block
03486 do i = 1, nx_block
03487 if (n == 2) sss (i,j,:) = c0
03488 if (n == 3) hmix (i,j,:) = c0
03489 if (n == 4) uocn (i,j,:) = c0
03490 if (n == 5) vocn (i,j,:) = c0
03491 if (n == 6) ss_tltx(i,j,:) = c0
03492 if (n == 7) ss_tlty(i,j,:) = c0
03493 if (n == 8) qdp (i,j,:) = c0
03494 do iblk = 1, nblocks
03495 if (hm(i,j,iblk) == c1) then
03496 if (n == 2) sss (i,j,iblk) = work1(i,j,iblk)
03497 if (n == 3) hmix (i,j,iblk) = work1(i,j,iblk)
03498 if (n == 4) uocn (i,j,iblk) = work1(i,j,iblk)
03499 if (n == 5) vocn (i,j,iblk) = work1(i,j,iblk)
03500 if (n == 6) ss_tltx(i,j,iblk) = work1(i,j,iblk)
03501 if (n == 7) ss_tlty(i,j,iblk) = work1(i,j,iblk)
03502 if (n == 8) qdp (i,j,iblk) = work1(i,j,iblk)
03503 endif
03504 enddo
03505 enddo
03506 enddo
03507 enddo
03508
03509 do j = 1, ny_block
03510 do i = 1, nx_block
03511 sss (i,j,:) = max (sss(i,j,:), c0)
03512 if (trim(Tfrzpt) == 'constant') then
03513 Tf (i,j,:) = -1.8_dbl_kind
03514 else
03515 Tf (i,j,:) = -depressT * sss(i,j,:)
03516 endif
03517 hmix(i,j,:) = max(hmix(i,j,:), c0)
03518 enddo
03519 enddo
03520
03521 if (restore_sst) then
03522 do j = 1, ny_block
03523 do i = 1, nx_block
03524 sst(i,j,:) = sst(i,j,:) + (work1(i,j,:)-sst(i,j,:))*dt/trest
03525 enddo
03526 enddo
03527
03528 endif
03529
03530
03531 if (istep1 <= 1 .and. .not. (restart)) then
03532 call interpolate_data (sst_data,sst)
03533 do iblk = 1, nblocks
03534 do j = 1, ny_block
03535 do i = 1, nx_block
03536 if (hm(i,j,iblk) == c1) then
03537 sst(i,j,iblk) = max (sst(i,j,iblk), Tf(i,j,iblk))
03538 else
03539 sst(i,j,iblk) = c0
03540 endif
03541 enddo
03542 enddo
03543 enddo
03544 endif
03545
03546 if (dbug) then
03547 if (my_task == master_task) &
03548 write (nu_diag,*) 'ocn_data_ncar'
03549 vmin = global_minval(Tf,distrb_info,tmask)
03550 vmax = global_maxval(Tf,distrb_info,tmask)
03551 if (my_task.eq.master_task) &
03552 write (nu_diag,*) 'Tf',vmin,vmax
03553 vmin = global_minval(sst,distrb_info,tmask)
03554 vmax = global_maxval(sst,distrb_info,tmask)
03555 if (my_task.eq.master_task) &
03556 write (nu_diag,*) 'sst',vmin,vmax
03557 vmin = global_minval(sss,distrb_info,tmask)
03558 vmax = global_maxval(sss,distrb_info,tmask)
03559 if (my_task.eq.master_task) &
03560 write (nu_diag,*) 'sss',vmin,vmax
03561 vmin = global_minval(hmix,distrb_info,tmask)
03562 vmax = global_maxval(hmix,distrb_info,tmask)
03563 if (my_task.eq.master_task) &
03564 write (nu_diag,*) 'hmix',vmin,vmax
03565 vmin = global_minval(uocn,distrb_info,umask)
03566 vmax = global_maxval(uocn,distrb_info,umask)
03567 if (my_task.eq.master_task) &
03568 write (nu_diag,*) 'uocn',vmin,vmax
03569 vmin = global_minval(vocn,distrb_info,umask)
03570 vmax = global_maxval(vocn,distrb_info,umask)
03571 if (my_task.eq.master_task) &
03572 write (nu_diag,*) 'vocn',vmin,vmax
03573 vmin = global_minval(ss_tltx,distrb_info,umask)
03574 vmax = global_maxval(ss_tltx,distrb_info,umask)
03575 if (my_task.eq.master_task) &
03576 write (nu_diag,*) 'ss_tltx',vmin,vmax
03577 vmin = global_minval(ss_tlty,distrb_info,umask)
03578 vmax = global_maxval(ss_tlty,distrb_info,umask)
03579 if (my_task.eq.master_task) &
03580 write (nu_diag,*) 'ss_tlty',vmin,vmax
03581 vmin = global_minval(qdp,distrb_info,tmask)
03582 vmax = global_maxval(qdp,distrb_info,tmask)
03583 if (my_task.eq.master_task) &
03584 write (nu_diag,*) 'qdp',vmin,vmax
03585 endif
03586
03587 end subroutine ocn_data_ncar
03588
03589
03590
03591
03592
03593
03594
03595
03596
03597 subroutine ocn_data_hadgem(dt)
03598
03599
03600
03601
03602
03603
03604
03605
03606
03607
03608
03609
03610
03611
03612 use ice_domain, only: nblocks
03613 use ice_flux, only: sst, uocn, vocn
03614 use ice_grid, only: t2ugrid_vector, ANGLET
03615
03616
03617
03618
03619
03620
03621 real (kind=dbl_kind), intent(in) ::
03622 dt
03623
03624 integer (kind=int_kind) ::
03625 i, j ,
03626 n ,
03627 iblk ,
03628 ixm,ixx,ixp ,
03629 recnum ,
03630 maxrec ,
03631 recslot ,
03632 dataloc ,
03633
03634 midmonth
03635
03636 real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) ::
03637 sstdat
03638
03639 real (kind=dbl_kind) :: workx, worky
03640
03641 logical (kind=log_kind) :: readm
03642
03643 character (char_len) ::
03644 fieldname
03645
03646 character (char_len_long) ::
03647 filename
03648
03649
03650
03651
03652
03653
03654
03655
03656 midmonth = 15
03657
03658
03659
03660 maxrec = 12
03661 ixm = mod(month+maxrec-2,maxrec) + 1
03662 ixp = mod(month, maxrec) + 1
03663 if (mday >= midmonth) ixm = 99
03664 if (mday < midmonth) ixp = 99
03665
03666
03667
03668
03669
03670
03671 recslot = 1
03672 if (mday < midmonth) recslot = 2
03673
03674
03675 call interp_coeff_monthly (recslot)
03676
03677
03678 readm = .false.
03679 if (istep==1 .or. (mday==midmonth .and. sec==0)) readm = .true.
03680
03681
03682 if (my_task == master_task .and. istep == 1) then
03683 write (nu_diag,*) ' '
03684 write (nu_diag,*) 'SST data interpolated to timestep:'
03685 write (nu_diag,*) trim(ocn_data_dir)//'MONTHLY/sst.1997.nc'
03686 if (restore_sst) write (nu_diag,*) &
03687 'SST restoring timescale (days) =', trestore
03688 if (trim(sst_data_type)=='hadgem_sst_uvocn') then
03689 write (nu_diag,*) ' '
03690 write (nu_diag,*) 'uocn and vocn interpolated to timestep:'
03691 write (nu_diag,*) trim(ocn_data_dir)//'MONTHLY/uocn.1997.nc'
03692 write (nu_diag,*) trim(ocn_data_dir)//'MONTHLY/vocn.1997.nc'
03693 endif
03694 endif
03695
03696
03697
03698
03699
03700 sst_file = trim(ocn_data_dir)//'MONTHLY/sst.1997.nc'
03701 fieldname='sst'
03702 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
03703 maxrec, sst_file, fieldname, sst_data, &
03704 field_loc_center, field_type_scalar)
03705
03706
03707 call interpolate_data (sst_data, sstdat)
03708
03709
03710 if (restore_sst) then
03711 do iblk = 1, nblocks
03712 do j = 1, ny_block
03713 do i = 1, nx_block
03714 sst(i,j,iblk) = sst(i,j,iblk) &
03715 + (sstdat(i,j,iblk)-sst(i,j,iblk))*dt/trest
03716 enddo
03717 enddo
03718 enddo
03719 endif
03720
03721
03722
03723
03724
03725
03726
03727
03728
03729
03730
03731 if (trim(sst_data_type)=='hadgem_sst_uvocn') then
03732
03733 filename = trim(ocn_data_dir)//'MONTHLY/uocn.1997.nc'
03734 fieldname='uocn'
03735 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
03736 maxrec, filename, fieldname, uocn_data, &
03737 field_loc_center, field_type_vector)
03738
03739
03740 call interpolate_data (uocn_data, uocn)
03741
03742 filename = trim(ocn_data_dir)//'MONTHLY/vocn.1997.nc'
03743 fieldname='vocn'
03744 call read_data_nc (readm, 0, fyear, ixm, month, ixp, &
03745 maxrec, filename, fieldname, vocn_data, &
03746 field_loc_center, field_type_vector)
03747
03748
03749 call interpolate_data (vocn_data, vocn)
03750
03751
03752
03753
03754
03755
03756 do iblk = 1, nblocks
03757 do j = 1, ny_block
03758 do i = 1, nx_block
03759
03760 workx = uocn(i,j,iblk)
03761 worky = vocn(i,j,iblk)
03762 uocn(i,j,iblk) = workx*cos(ANGLET(i,j,iblk)) &
03763 + worky*sin(ANGLET(i,j,iblk))
03764 vocn(i,j,iblk) = worky*cos(ANGLET(i,j,iblk)) &
03765 - workx*sin(ANGLET(i,j,iblk))
03766
03767 uocn(i,j,iblk) = uocn(i,j,iblk) * cm_to_m
03768 vocn(i,j,iblk) = vocn(i,j,iblk) * cm_to_m
03769
03770 enddo
03771 enddo
03772 enddo
03773
03774
03775
03776
03777
03778 call t2ugrid_vector(uocn)
03779 call t2ugrid_vector(vocn)
03780
03781
03782 endif
03783
03784
03785 end subroutine ocn_data_hadgem
03786
03787
03788
03789
03790 end module ice_forcing
03791
03792