Newer
Older
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn
sdf(jf)%igrd = 0
sdf(jf)%ibdy = 0
sdf(jf)%imap => NULL()
sdf(jf)%ltotvel = .FALSE.
sdf(jf)%lzint = .FALSE.
END DO
!
IF(lwp) THEN ! control print
WRITE(numout,*)
IF( .NOT.PRESENT( knoprint) ) THEN
WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
ENDIF
WRITE(numout,*) ' fld_fill : fill data structure with information from namelist '//TRIM( cdnam )
WRITE(numout,*) ' ~~~~~~~~'
WRITE(numout,*) ' list of files and frequency (>0: in hours ; <0 in months)'
DO jf = 1, SIZE(sdf)
WRITE(numout,*) ' root filename: ' , TRIM( sdf(jf)%clrootname ), ' variable name: ', TRIM( sdf(jf)%clvar )
WRITE(numout,*) ' frequency: ' , sdf(jf)%freqh , &
& ' time interp: ' , sdf(jf)%ln_tint , &
& ' climatology: ' , sdf(jf)%ln_clim
WRITE(numout,*) ' weights: ' , TRIM( sdf(jf)%wgtname ), &
& ' pairing: ' , TRIM( sdf(jf)%vcomp ), &
& ' data type: ' , sdf(jf)%clftyp , &
& ' land/sea mask:' , TRIM( sdf(jf)%lsmname )
call flush(numout)
END DO
ENDIF
!
END SUBROUTINE fld_fill
SUBROUTINE wgt_list( sd, kwgt )
!!---------------------------------------------------------------------
!! *** ROUTINE wgt_list ***
!!
!! ** Purpose : search array of WGTs and find a weights file entry,
!! or return a new one adding it to the end if new entry.
!! the weights data is read in and restructured (fld_weight)
!!----------------------------------------------------------------------
TYPE( FLD ), INTENT(in ) :: sd ! field with name of weights file
INTEGER , INTENT( out) :: kwgt ! index of weights
!
INTEGER :: kw, nestid ! local integer
!!----------------------------------------------------------------------
!
!! search down linked list
!! weights filename is either present or we hit the end of the list
!
!! because agrif nest part of filenames are now added in iom_open
!! to distinguish between weights files on the different grids, need to track
!! nest number explicitly
nestid = 0
#if defined key_agrif
nestid = Agrif_Fixed()
#endif
DO kw = 1, nxt_wgt-1
IF( ref_wgts(kw)%wgtname == sd%wgtname .AND. &
ref_wgts(kw)%nestid == nestid) THEN
kwgt = kw
RETURN
ENDIF
END DO
kwgt = nxt_wgt
CALL fld_weight( sd )
!
END SUBROUTINE wgt_list
SUBROUTINE wgt_print( )
!!---------------------------------------------------------------------
!! *** ROUTINE wgt_print ***
!!
!! ** Purpose : print the list of known weights
!!----------------------------------------------------------------------
INTEGER :: kw !
!!----------------------------------------------------------------------
!
DO kw = 1, nxt_wgt-1
WRITE(numout,*) 'weight file: ',TRIM(ref_wgts(kw)%wgtname)
WRITE(numout,*) ' ddims: ',ref_wgts(kw)%ddims(1),ref_wgts(kw)%ddims(2)
WRITE(numout,*) ' numwgt: ',ref_wgts(kw)%numwgt
WRITE(numout,*) ' jpiwgt: ',ref_wgts(kw)%jpiwgt
WRITE(numout,*) ' jpjwgt: ',ref_wgts(kw)%jpjwgt
WRITE(numout,*) ' botleft: ',ref_wgts(kw)%botleft
WRITE(numout,*) ' topright: ',ref_wgts(kw)%topright
IF( ref_wgts(kw)%cyclic ) THEN
WRITE(numout,*) ' cyclical'
IF( ref_wgts(kw)%overlap > 0 ) WRITE(numout,*) ' with overlap of ', ref_wgts(kw)%overlap
ELSE
WRITE(numout,*) ' not cyclical'
ENDIF
IF( ASSOCIATED(ref_wgts(kw)%data_wgt) ) WRITE(numout,*) ' allocated'
END DO
!
END SUBROUTINE wgt_print
SUBROUTINE fld_weight( sd )
!!---------------------------------------------------------------------
!! *** ROUTINE fld_weight ***
!!
!! ** Purpose : create a new WGT structure and fill in data from file,
!! restructuring as required
!!----------------------------------------------------------------------
TYPE( FLD ), INTENT(in) :: sd ! field with name of weights file
!!
INTEGER :: ji,jj,jn ! dummy loop indices
INTEGER :: inum ! local logical unit
INTEGER :: id ! local variable id
INTEGER :: ipk ! local vertical dimension
INTEGER :: zwrap ! local integer
LOGICAL :: cyclical !
CHARACTER (len=5) :: clname !
INTEGER , DIMENSION(4) :: ddims
INTEGER :: isrc
sparonuz
committed
REAL(dp), DIMENSION(jpi,jpj) :: data_tmp
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
!!----------------------------------------------------------------------
!
IF( nxt_wgt > tot_wgts ) THEN
CALL ctl_stop("fld_weight: weights array size exceeded, increase tot_wgts")
ENDIF
!
!! new weights file entry, add in extra information
!! a weights file represents a 2D grid of a certain shape, so we assume that the current
!! input data file is representative of all other files to be opened and processed with the
!! current weights file
!! get data grid dimensions
id = iom_varid( sd%num, sd%clvar, ddims )
!! now open the weights file
CALL iom_open ( sd%wgtname, inum ) ! interpolation weights
IF( inum > 0 ) THEN
!! determine whether we have an east-west cyclic grid
!! from global attribute called "ew_wrap" in the weights file
!! note that if not found, iom_getatt returns -999 and cyclic with no overlap is assumed
!! since this is the most common forcing configuration
CALL iom_getatt(inum, 'ew_wrap', zwrap)
IF( zwrap >= 0 ) THEN
cyclical = .TRUE.
ELSE IF( zwrap == -999 ) THEN
cyclical = .TRUE.
zwrap = 0
ELSE
cyclical = .FALSE.
ENDIF
ref_wgts(nxt_wgt)%ddims(1) = ddims(1)
ref_wgts(nxt_wgt)%ddims(2) = ddims(2)
ref_wgts(nxt_wgt)%wgtname = sd%wgtname
ref_wgts(nxt_wgt)%overlap = zwrap
ref_wgts(nxt_wgt)%cyclic = cyclical
ref_wgts(nxt_wgt)%nestid = 0
#if defined key_agrif
ref_wgts(nxt_wgt)%nestid = Agrif_Fixed()
#endif
!! weights file is stored as a set of weights (wgt01->wgt04 or wgt01->wgt16)
!! for each weight wgtNN there is an integer array srcNN which gives the point in
!! the input data grid which is to be multiplied by the weight
!! they are both arrays on the model grid so the result of the multiplication is
!! added into an output array on the model grid as a running sum
!! two possible cases: bilinear (4 weights) or bicubic (16 weights)
id = iom_varid(inum, 'src05', ldstop=.FALSE.)
IF( id <= 0 ) THEN ; ref_wgts(nxt_wgt)%numwgt = 4
ELSE ; ref_wgts(nxt_wgt)%numwgt = 16
ENDIF
ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(Nis0:Nie0,Njs0:Nje0,4) )
ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(Nis0:Nie0,Njs0:Nje0,4) )
ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(Nis0:Nie0,Njs0:Nje0,ref_wgts(nxt_wgt)%numwgt) )
DO jn = 1,4
WRITE(clname,'(a3,i2.2)') 'src',jn
CALL iom_get ( inum, jpdom_global, clname, data_tmp(:,:), cd_type = 'Z' ) ! no call to lbc_lnk
DO_2D( 0, 0, 0, 0 )
isrc = NINT(data_tmp(ji,jj)) - 1
ref_wgts(nxt_wgt)%data_jpi(ji,jj,jn) = 1 + MOD(isrc, ref_wgts(nxt_wgt)%ddims(1))
ref_wgts(nxt_wgt)%data_jpj(ji,jj,jn) = 1 + isrc / ref_wgts(nxt_wgt)%ddims(1)
END_2D
END DO
DO jn = 1, ref_wgts(nxt_wgt)%numwgt
WRITE(clname,'(a3,i2.2)') 'wgt',jn
CALL iom_get ( inum, jpdom_global, clname, data_tmp(:,:), cd_type = 'Z' ) ! no call to lbc_lnk
DO_2D( 0, 0, 0, 0 )
ref_wgts(nxt_wgt)%data_wgt(ji,jj,jn) = data_tmp(ji,jj)
END_2D
END DO
CALL iom_close (inum)
! find min and max indices in grid
ref_wgts(nxt_wgt)%botleft( 1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
ref_wgts(nxt_wgt)%botleft( 2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:))
ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:))
! and therefore dimensions of the input box
ref_wgts(nxt_wgt)%jpiwgt = ref_wgts(nxt_wgt)%topright(1) - ref_wgts(nxt_wgt)%botleft(1) + 1
ref_wgts(nxt_wgt)%jpjwgt = ref_wgts(nxt_wgt)%topright(2) - ref_wgts(nxt_wgt)%botleft(2) + 1
! shift indexing of source grid
ref_wgts(nxt_wgt)%data_jpi(:,:,:) = ref_wgts(nxt_wgt)%data_jpi(:,:,:) - ref_wgts(nxt_wgt)%botleft(1) + 1
ref_wgts(nxt_wgt)%data_jpj(:,:,:) = ref_wgts(nxt_wgt)%data_jpj(:,:,:) - ref_wgts(nxt_wgt)%botleft(2) + 1
! create input grid, give it a halo to allow gradient calculations
! SA: +3 stencil is a patch to avoid out-of-bound computation in some configuration.
! a more robust solution will be given in next release
ipk = SIZE(sd%fnow, 3)
ALLOCATE( ref_wgts(nxt_wgt)%fly_dta(ref_wgts(nxt_wgt)%jpiwgt+3, ref_wgts(nxt_wgt)%jpjwgt+3 ,ipk) )
IF( ref_wgts(nxt_wgt)%cyclic ) ALLOCATE( ref_wgts(nxt_wgt)%col(1,ref_wgts(nxt_wgt)%jpjwgt+3,ipk) )
!
nxt_wgt = nxt_wgt + 1
!
ELSE
CALL ctl_stop( ' fld_weight : unable to read the file ' )
ENDIF
!
END SUBROUTINE fld_weight
SUBROUTINE apply_seaoverland( clmaskfile, zfieldo, jpi1_lsm, jpi2_lsm, jpj1_lsm, &
& jpj2_lsm, itmpi, itmpj, itmpz, rec1_lsm, recn_lsm )
!!---------------------------------------------------------------------
!! *** ROUTINE apply_seaoverland ***
!!
!! ** Purpose : avoid spurious fluxes in coastal or near-coastal areas
!! due to the wrong usage of "land" values from the coarse
!! atmospheric model when spatial interpolation is required
!! D. Delrosso INGV
!!----------------------------------------------------------------------
INTEGER, INTENT(in ) :: itmpi,itmpj,itmpz ! lengths
INTEGER, INTENT(in ) :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices
INTEGER, DIMENSION(3), INTENT(in ) :: rec1_lsm,recn_lsm ! temporary arrays for start and length
REAL(wp),DIMENSION (:,:,:),INTENT(inout) :: zfieldo ! input/output array for seaoverland application
CHARACTER (len=100), INTENT(in ) :: clmaskfile ! land/sea mask file name
!
INTEGER :: inum,jni,jnj,jnz,jc ! local indices
REAL(wp),DIMENSION (:,:,:),ALLOCATABLE :: zslmec1 ! local array for land point detection
REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfieldn ! array of forcing field with undeff for land points
REAL(wp),DIMENSION (:,:), ALLOCATABLE :: zfield ! array of forcing field
!!---------------------------------------------------------------------
!
ALLOCATE ( zslmec1(itmpi,itmpj,itmpz), zfieldn(itmpi,itmpj), zfield(itmpi,itmpj) )
!
! Retrieve the land sea mask data
CALL iom_open( clmaskfile, inum )
SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
CASE(1)
CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), &
& 1, kstart = rec1_lsm, kcount = recn_lsm)
CASE DEFAULT
CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), &
& 1, kstart = rec1_lsm, kcount = recn_lsm)
END SELECT
CALL iom_close( inum )
!
DO jnz=1,rec1_lsm(3) !! Loop over k dimension
!
DO jni = 1, itmpi !! copy the original field into a tmp array
DO jnj = 1, itmpj !! substituting undeff over land points
zfieldn(jni,jnj) = zfieldo(jni,jnj,jnz)
IF( zslmec1(jni,jnj,jnz) == 1. ) zfieldn(jni,jnj) = undeff_lsm
END DO
END DO
!
CALL seaoverland( zfieldn, itmpi, itmpj, zfield )
DO jc = 1, nn_lsm
CALL seaoverland( zfield, itmpi, itmpj, zfield )
END DO
!
! Check for Undeff and substitute original values
IF( ANY(zfield==undeff_lsm) ) THEN
DO jni = 1, itmpi
DO jnj = 1, itmpj
IF( zfield(jni,jnj)==undeff_lsm ) zfield(jni,jnj) = zfieldo(jni,jnj,jnz)
END DO
END DO
ENDIF
!
zfieldo(:,:,jnz) = zfield(:,:)
!
END DO !! End Loop over k dimension
!
DEALLOCATE ( zslmec1, zfieldn, zfield )
!
END SUBROUTINE apply_seaoverland
SUBROUTINE seaoverland( zfieldn, ileni, ilenj, zfield )
!!---------------------------------------------------------------------
!! *** ROUTINE seaoverland ***
!!
!! ** Purpose : create shifted matrices for seaoverland application
!! D. Delrosso INGV
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: ileni,ilenj ! lengths
REAL(wp), DIMENSION (ileni,ilenj), INTENT(in ) :: zfieldn ! array of forcing field with undeff for land points
REAL(wp), DIMENSION (ileni,ilenj), INTENT( out) :: zfield ! array of forcing field
!
REAL(wp) , DIMENSION (ileni,ilenj) :: zmat1, zmat2, zmat3, zmat4 ! local arrays
REAL(wp) , DIMENSION (ileni,ilenj) :: zmat5, zmat6, zmat7, zmat8 ! - -
REAL(wp) , DIMENSION (ileni,ilenj) :: zlsm2d ! - -
REAL(wp) , DIMENSION (ileni,ilenj,8) :: zlsm3d ! - -
LOGICAL , DIMENSION (ileni,ilenj,8) :: ll_msknan3d ! logical mask for undeff detection
LOGICAL , DIMENSION (ileni,ilenj) :: ll_msknan2d ! logical mask for undeff detection
!!----------------------------------------------------------------------
zmat8 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(:,1)/) , DIM=2 )
zmat1 = eoshift( zmat8 , SHIFT=-1 , BOUNDARY = (/zmat8(1,:)/) , DIM=1 )
zmat2 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(1,:)/) , DIM=1 )
zmat4 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(:,ilenj)/) , DIM=2 )
zmat3 = eoshift( zmat4 , SHIFT=-1 , BOUNDARY = (/zmat4(1,:)/) , DIM=1 )
zmat5 = eoshift( zmat4 , SHIFT= 1 , BOUNDARY = (/zmat4(ileni,:)/) , DIM=1 )
zmat6 = eoshift( zfieldn , SHIFT= 1 , BOUNDARY = (/zfieldn(ileni,:)/) , DIM=1 )
zmat7 = eoshift( zmat8 , SHIFT= 1 , BOUNDARY = (/zmat8(ileni,:)/) , DIM=1 )
!
zlsm3d = RESHAPE( (/ zmat1, zmat2, zmat3, zmat4, zmat5, zmat6, zmat7, zmat8 /), (/ ileni, ilenj, 8 /))
ll_msknan3d = .NOT.( zlsm3d == undeff_lsm )
ll_msknan2d = .NOT.( zfieldn == undeff_lsm ) ! FALSE where is Undeff (land)
zlsm2d = SUM( zlsm3d, 3 , ll_msknan3d ) / MAX( 1 , COUNT( ll_msknan3d , 3 ) )
WHERE( COUNT( ll_msknan3d , 3 ) == 0._wp ) zlsm2d = undeff_lsm
zfield = MERGE( zfieldn, zlsm2d, ll_msknan2d )
!
END SUBROUTINE seaoverland
SUBROUTINE fld_interp( num, clvar, kw, kk, dta, nrec, lsmfile)
!!---------------------------------------------------------------------
!! *** ROUTINE fld_interp ***
!!
!! ** Purpose : apply weights to input gridded data to create data
!! on model grid
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: num ! stream number
CHARACTER(LEN=*) , INTENT(in ) :: clvar ! variable name
INTEGER , INTENT(in ) :: kw ! weights number
INTEGER , INTENT(in ) :: kk ! vertical dimension of kk
REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: dta ! output field on model grid
INTEGER , INTENT(in ) :: nrec ! record number to read (ie time slice)
CHARACTER(LEN=*) , INTENT(in ) :: lsmfile ! land sea mask file name
!
INTEGER, DIMENSION(3) :: rec1, recn ! temporary arrays for start and length
INTEGER, DIMENSION(3) :: rec1_lsm, recn_lsm ! temporary arrays for start and length in case of seaoverland
INTEGER :: ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2 ! temporary indices
INTEGER :: ji, jj, jk, jn, jir, jjr ! loop counters
INTEGER :: ipk
INTEGER :: ni, nj ! lengths
INTEGER :: jpimin,jpiwid ! temporary indices
INTEGER :: jpimin_lsm,jpiwid_lsm ! temporary indices
INTEGER :: jpjmin,jpjwid ! temporary indices
INTEGER :: jpjmin_lsm,jpjwid_lsm ! temporary indices
INTEGER :: jpi1,jpi2,jpj1,jpj2 ! temporary indices
INTEGER :: jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm ! temporary indices
INTEGER :: itmpi,itmpj,itmpz ! lengths
REAL(wp),DIMENSION(:,:,:), ALLOCATABLE :: ztmp_fly_dta ! local array of values on input grid
!!----------------------------------------------------------------------
ipk = SIZE(dta, 3)
!
!! for weighted interpolation we have weights at four corners of a box surrounding
!! a model grid point, each weight is multiplied by a grid value (bilinear case)
!! or by a grid value and gradients at the corner point (bicubic case)
!! so we need to have a 4 by 4 subgrid surrounding each model point to cover both cases
!! sub grid from non-model input grid which encloses all grid points in this nemo process
jpimin = ref_wgts(kw)%botleft(1)
jpjmin = ref_wgts(kw)%botleft(2)
jpiwid = ref_wgts(kw)%jpiwgt
jpjwid = ref_wgts(kw)%jpjwgt
!! when reading in, expand this sub-grid by one halo point all the way round for calculating gradients
rec1(1) = MAX( jpimin-1, 1 )
rec1(2) = MAX( jpjmin-1, 1 )
rec1(3) = 1
recn(1) = MIN( jpiwid+2, ref_wgts(kw)%ddims(1)-rec1(1)+1 )
recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
recn(3) = kk
!! where we need to put it in the non-nemo grid fly_dta
!! note that jpi1 and jpj1 only differ from 1 when jpimin and jpjmin are 1
!! (ie at the extreme west or south of the whole input grid) and similarly for jpi2 and jpj2
jpi1 = 2 + rec1(1) - jpimin
jpj1 = 2 + rec1(2) - jpjmin
jpi2 = jpi1 + recn(1) - 1
jpj2 = jpj1 + recn(2) - 1
IF( LEN_TRIM(lsmfile) > 0 ) THEN
!! indeces for ztmp_fly_dta
! --------------------------
rec1_lsm(1)=MAX(rec1(1)-nn_lsm,1) ! starting index for enlarged external data, x direction
rec1_lsm(2)=MAX(rec1(2)-nn_lsm,1) ! starting index for enlarged external data, y direction
rec1_lsm(3) = 1 ! vertical dimension
recn_lsm(1)=MIN(rec1(1)-rec1_lsm(1)+recn(1)+nn_lsm,ref_wgts(kw)%ddims(1)-rec1_lsm(1)) ! n points in x direction
recn_lsm(2)=MIN(rec1(2)-rec1_lsm(2)+recn(2)+nn_lsm,ref_wgts(kw)%ddims(2)-rec1_lsm(2)) ! n points in y direction
recn_lsm(3) = kk ! number of vertical levels in the input file
! Avoid out of bound
jpimin_lsm = MAX( rec1_lsm(1)+1, 1 )
jpjmin_lsm = MAX( rec1_lsm(2)+1, 1 )
jpiwid_lsm = MIN( recn_lsm(1)-2,ref_wgts(kw)%ddims(1)-rec1(1)+1)
jpjwid_lsm = MIN( recn_lsm(2)-2,ref_wgts(kw)%ddims(2)-rec1(2)+1)
jpi1_lsm = 2+rec1_lsm(1)-jpimin_lsm
jpj1_lsm = 2+rec1_lsm(2)-jpjmin_lsm
jpi2_lsm = jpi1_lsm + recn_lsm(1) - 1
jpj2_lsm = jpj1_lsm + recn_lsm(2) - 1
itmpi=jpi2_lsm-jpi1_lsm+1
itmpj=jpj2_lsm-jpj1_lsm+1
itmpz=kk
ALLOCATE(ztmp_fly_dta(itmpi,itmpj,itmpz))
ztmp_fly_dta(:,:,:) = 0.0
SELECT CASE( SIZE(ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) )
CASE(1)
CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), &
& nrec, kstart = rec1_lsm, kcount = recn_lsm)
CASE DEFAULT
CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), &
& nrec, kstart = rec1_lsm, kcount = recn_lsm)
END SELECT
CALL apply_seaoverland(lsmfile,ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), &
& jpi1_lsm,jpi2_lsm,jpj1_lsm,jpj2_lsm, &
& itmpi,itmpj,itmpz,rec1_lsm,recn_lsm)
! Relative indeces for remapping
ii_lsm1 = (rec1(1)-rec1_lsm(1))+1
ii_lsm2 = (ii_lsm1+recn(1))-1
ij_lsm1 = (rec1(2)-rec1_lsm(2))+1
ij_lsm2 = (ij_lsm1+recn(2))-1
ref_wgts(kw)%fly_dta(:,:,:) = 0.0
ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:) = ztmp_fly_dta(ii_lsm1:ii_lsm2,ij_lsm1:ij_lsm2,:)
DEALLOCATE(ztmp_fly_dta)
ELSE
ref_wgts(kw)%fly_dta(:,:,:) = 0.0
CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, kstart = rec1, kcount = recn)
ENDIF
!! first four weights common to both bilinear and bicubic
!! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft
!! note that we have to offset by 1 into fly_dta array because of halo added to fly_dta (rec1 definition)
dta(:,:,:) = 0._wp
DO jn = 1,4
DO_3D( 0, 0, 0, 0, 1,ipk )
ni = ref_wgts(kw)%data_jpi(ji,jj,jn) + 1
nj = ref_wgts(kw)%data_jpj(ji,jj,jn) + 1
dta(ji,jj,jk) = dta(ji,jj,jk) + ref_wgts(kw)%data_wgt(ji,jj,jn) * ref_wgts(kw)%fly_dta(ni,nj,jk)
END_3D
END DO
IF(ref_wgts(kw)%numwgt .EQ. 16) THEN
!! fix up halo points that we couldnt read from file
IF( jpi1 == 2 ) THEN
ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:)
ENDIF
IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:)
ENDIF
IF( jpj1 == 2 ) THEN
ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:)
ENDIF
IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .LT. jpjwid+2 ) THEN
ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:)
ENDIF
!! if data grid is cyclic we can do better on east-west edges
!! but have to allow for whether first and last columns are coincident
IF( ref_wgts(kw)%cyclic ) THEN
rec1(2) = MAX( jpjmin-1, 1 )
recn(1) = 1
recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 )
jpj1 = 2 + rec1(2) - jpjmin
jpj2 = jpj1 + recn(2) - 1
IF( jpi1 == 2 ) THEN
rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap
CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, kstart = rec1, kcount = recn)
ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
ENDIF
IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN
rec1(1) = 1 + ref_wgts(kw)%overlap
CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, kstart = rec1, kcount = recn)
ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:)
ENDIF
ENDIF
!
!!$ DO jn = 1,4
!!$ DO_3D( 0, 0, 0, 0, 1,ipk )
!!$ ni = ref_wgts(kw)%data_jpi(ji,jj,jn) + 1
!!$ nj = ref_wgts(kw)%data_jpj(ji,jj,jn) + 1
!!$ dta(ji,jj,jk) = dta(ji,jj,jk) &
!!$ ! gradient in the i direction
!!$ & + ref_wgts(kw)%data_wgt(ji,jj,jn+4) * 0.5_wp * &
!!$ & (ref_wgts(kw)%fly_dta(ni+1,nj ,jk) - ref_wgts(kw)%fly_dta(ni-1,nj ,jk)) &
!!$ ! gradient in the j direction
!!$ & + ref_wgts(kw)%data_wgt(ji,jj,jn+8) * 0.5_wp * &
!!$ & (ref_wgts(kw)%fly_dta(ni ,nj+1,jk) - ref_wgts(kw)%fly_dta(ni ,nj-1,jk)) &
!!$ ! gradient in the ij direction
!!$ & + ref_wgts(kw)%data_wgt(ji,jj,jn+12) * 0.25_wp * &
!!$ & ((ref_wgts(kw)%fly_dta(ni+1,nj+1,jk) - ref_wgts(kw)%fly_dta(ni-1,nj+1,jk)) - &
!!$ & (ref_wgts(kw)%fly_dta(ni+1,nj-1,jk) - ref_wgts(kw)%fly_dta(ni-1,nj-1,jk)))
!!$ END_3D
!!$ END DO
!
DO jn = 1,4
DO_3D( 0, 0, 0, 0, 1,ipk )
ni = ref_wgts(kw)%data_jpi(ji,jj,jn)
nj = ref_wgts(kw)%data_jpj(ji,jj,jn)
! gradient in the i direction
dta(ji,jj,jk) = dta(ji,jj,jk) + ref_wgts(kw)%data_wgt(ji,jj,jn+4) * 0.5_wp * &
& (ref_wgts(kw)%fly_dta(ni+2,nj+1,jk) - ref_wgts(kw)%fly_dta(ni ,nj+1,jk))
END_3D
END DO
DO jn = 1,4
DO_3D( 0, 0, 0, 0, 1,ipk )
ni = ref_wgts(kw)%data_jpi(ji,jj,jn)
nj = ref_wgts(kw)%data_jpj(ji,jj,jn)
! gradient in the j direction
dta(ji,jj,jk) = dta(ji,jj,jk) + ref_wgts(kw)%data_wgt(ji,jj,jn+8) * 0.5_wp * &
& (ref_wgts(kw)%fly_dta(ni+1,nj+2,jk) - ref_wgts(kw)%fly_dta(ni+1,nj ,jk))
END_3D
END DO
DO jn = 1,4
DO_3D( 0, 0, 0, 0, 1,ipk )
ni = ref_wgts(kw)%data_jpi(ji,jj,jn)
nj = ref_wgts(kw)%data_jpj(ji,jj,jn)
! gradient in the ij direction
dta(ji,jj,jk) = dta(ji,jj,jk) + ref_wgts(kw)%data_wgt(ji,jj,jn+12) * 0.25_wp * ( &
& (ref_wgts(kw)%fly_dta(ni+2,nj+2,jk) - ref_wgts(kw)%fly_dta(ni ,nj+2,jk)) - &
& (ref_wgts(kw)%fly_dta(ni+2,nj ,jk) - ref_wgts(kw)%fly_dta(ni ,nj ,jk)))
END_3D
END DO
!
ENDIF
!
END SUBROUTINE fld_interp
FUNCTION fld_filename( sdjf, kday, kmonth, kyear )
!!---------------------------------------------------------------------
!! *** FUNCTION fld_filename ***
!!
!! ** Purpose : define the filename according to a given date
!!---------------------------------------------------------------------
TYPE(FLD), INTENT(in) :: sdjf ! input field related variables
INTEGER , INTENT(in) :: kday, kmonth, kyear
!
CHARACTER(len = 256) :: clname, fld_filename
!!---------------------------------------------------------------------
! build the new filename if not climatological data
clname=TRIM(sdjf%clrootname)
!
! note that sdjf%ln_clim is is only acting on the presence of the year in the file name
IF( .NOT. sdjf%ln_clim ) THEN
WRITE(clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), kyear ! add year
IF( sdjf%clftyp /= 'yearly' ) WRITE(clname, '(a, "m",i2.2)' ) TRIM( clname ), kmonth ! add month
ELSE
! build the new filename if climatological data
IF( sdjf%clftyp /= 'yearly' ) WRITE(clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), kmonth ! add month
ENDIF
IF( sdjf%clftyp == 'daily' .OR. sdjf%clftyp(1:4) == 'week' ) &
& WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname ), kday ! add day
fld_filename = clname
END FUNCTION fld_filename
FUNCTION ksec_week( cdday )
!!---------------------------------------------------------------------
!! *** FUNCTION ksec_week ***
!!
!! ** Purpose : seconds between 00h of the beginning of the week and half of the current time step
!!---------------------------------------------------------------------
CHARACTER(len=*), INTENT(in) :: cdday ! first 3 letters of the first day of the weekly file
!!
INTEGER :: ksec_week ! output variable
INTEGER :: ijul, ishift ! local integer
CHARACTER(len=3),DIMENSION(7) :: cl_week
!!----------------------------------------------------------------------
cl_week = (/"sun","sat","fri","thu","wed","tue","mon"/)
DO ijul = 1, 7
IF( cl_week(ijul) == TRIM(cdday) ) EXIT
END DO
IF( ijul .GT. 7 ) CALL ctl_stop( 'ksec_week: wrong day for sdjf%clftyp(6:8): '//TRIM(cdday) )
!
ishift = ijul * NINT(rday)
!
ksec_week = nsec_monday + ishift
ksec_week = MOD( ksec_week, 7*NINT(rday) )
!
END FUNCTION ksec_week
!!======================================================================
END MODULE fldread