- Timestamp:
- 2022-11-08T16:24:46+01:00 (20 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/forcing_tools.f90
r7475 r7796 51 51 USE mod_orchidee_para 52 52 USE forcingdaily_tools 53 USE grid 53 54 ! 54 55 IMPLICIT NONE … … 90 91 REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: time_bounds 91 92 CHARACTER(LEN=20), SAVE, ALLOCATABLE, DIMENSION(:) :: time_axename, time_cellmethod 92 REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: preciptime93 93 INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: time_sourcefile 94 94 INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: time_id … … 404 404 ! 405 405 SUBROUTINE forcing_open(filenames_in, iim, jjm, lon, lat, nbpoint_in, drvzoom_lon, drvzoom_lat, & 406 & kindex, nbindex_perproc, wunit, landonly_arg)406 & kindex, nbindex_perproc, wunit, model_guess, landonly_arg) 407 407 ! 408 408 ! Opens the forcing file and reads some key information and stores them in the shared variables of the … … 424 424 INTEGER(i_std), OPTIONAL, INTENT(in) :: wunit 425 425 LOGICAL, OPTIONAL, INTENT(in) :: landonly_arg 426 CHARACTER(LEN=*), INTENT(in) :: model_guess 426 427 ! 427 428 ! LOCAL … … 501 502 ! 502 503 IF ( PRESENT(wunit) ) THEN 503 WRITE(wunit,*) "Getting the zoomedgrid", nbpoint_tmp504 WRITE(wunit,*) "Getting the full grid", nbpoint_tmp 504 505 CALL FLUSH(wunit) 505 506 ENDIF 506 CALL forcing_zoomgrid(drvzoom_lon, drvzoom_lat, forfilename(1), .FALSE.)507 CALL forcing_zoomgrid(drvzoom_lon, drvzoom_lat, forfilename(1), model_guess, .FALSE.) 507 508 IF ( PRESENT(wunit) ) THEN 508 509 WRITE(wunit,*) "Out of the zoomed grid operation" … … 1428 1429 ! Local 1429 1430 ! 1430 INTEGER(i_std) :: is 1431 INTEGER(i_std) :: is, tlen 1431 1432 REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tair_full, qair_full 1432 1433 REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: tairmin_full, tairmax_full … … 1435 1436 REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: u_full, v_full 1436 1437 REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: ps_full, ztq_full, zuv_full 1437 REAL(r_std), ALLOCATABLE, DIMENSION(:) :: preciptime_tmp 1438 REAL(r_std), ALLOCATABLE, DIMENSION(:) :: time_precip_old 1439 REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: savepreciptime 1440 INTEGER(i_std) :: s2in1(1) 1441 LOGICAL, SAVE :: first_call_readslab=.TRUE. 1438 1442 ! 1439 1443 ! 1.0 Verify that for the slabs the memory is allocated for the variable … … 1457 1461 IF ( .NOT. ALLOCATED(timebnd_precip) ) ALLOCATE(timebnd_precip(slab_size,2)) 1458 1462 IF ( .NOT. ALLOCATED(preciptime_slab) ) ALLOCATE(preciptime_slab(nbpoint_proc,slab_size)) 1459 IF ( .NOT. ALLOCATED(preciptime_tmp) ) ALLOCATE(preciptime_tmp(slab_size)) 1463 IF ( .NOT. ALLOCATED(time_precip_old) ) ALLOCATE(time_precip_old(slab_size)) 1464 IF ( .NOT. ALLOCATED(savepreciptime) ) ALLOCATE(savepreciptime(nbpoint_proc,slab_size)) 1465 1460 1466 ! 1461 1467 IF ( .NOT. ALLOCATED(swdown_slab) ) ALLOCATE(swdown_slab(nbpoint_proc,slab_size)) … … 1482 1488 IF ( .NOT. ALLOCATED(zuv_slab) ) ALLOCATE(zuv_slab(nbpoint_proc,slab_size)) 1483 1489 ! 1490 ! 1491 ! Save the old time axis for precipitation. 1492 ! 1493 IF ( first_call_readslab ) THEN 1494 time_precip_old(:) = zero 1495 preciptime_slab(:,:) = zero 1496 ELSE 1497 time_precip_old(:) = time_precip(:) 1498 ENDIF 1484 1499 ! 1485 1500 IF ( is_root_prc) THEN … … 1506 1521 ALLOCATE(zuv_full(nbpoint_loc,slab_size)) 1507 1522 ! 1508 CALL forcing_readslab_root(time_int, &1523 CALL forcing_readslab_root(time_int, first_call_readslab, & 1509 1524 & tair_full, tairmax_full, tairmin_full, time_tair, timebnd_tair, & 1510 1525 & qair_full, time_qair, timebnd_qair, & 1511 & rainf_full, snowf_full, time_precip, timebnd_precip, preciptime_tmp,&1526 & rainf_full, snowf_full, time_precip, timebnd_precip, & 1512 1527 & swdown_full, time_swdown, timebnd_swdown, & 1513 1528 & lwdown_full, time_lwdown, timebnd_lwdown, & … … 1546 1561 CALL bcast(time_precip) 1547 1562 CALL bcast(timebnd_precip) 1548 CALL bcast(preciptime_tmp)1549 1563 CALL bcast(time_swdown) 1550 1564 CALL bcast(timebnd_swdown) … … 1596 1610 ! 1597 1611 ! 1598 DO is=1,nbpoint_proc 1599 preciptime_slab(is,:) = preciptime_tmp(:) 1600 ENDDO 1612 ! Transfer preciptime_slab from the old slab to the new one. 1613 ! s2in1: the index of the starting time from new slab in the previous slab 1614 s2in1 = MINLOC(ABS(time_precip_old-time_precip(1))) 1615 IF ( size(s2in1) .EQ. 0 ) THEN 1616 s2in1(1) = 1 1617 ENDIF 1618 tlen=slab_size-s2in1(1)+1 1619 savepreciptime(:,:) = zero 1620 savepreciptime(:,1:tlen) = preciptime_slab(:,s2in1(1):slab_size) 1621 preciptime_slab(:,:) = zero 1622 preciptime_slab(:,1:tlen) = savepreciptime(:,1:tlen) 1623 1624 1625 IF (first_call_readslab) THEN 1626 first_call_readslab = .FALSE. 1627 ENDIF 1628 1601 1629 ! 1602 1630 ! Clean-up to free the memory on the root processor. … … 1630 1658 !_ ============================================================================================================================== 1631 1659 1632 SUBROUTINE forcing_readslab_root(time_int, &1660 SUBROUTINE forcing_readslab_root(time_int, first_call_readslab, & 1633 1661 & tair, tairmax, tairmin, t_tair, tbnd_tair, & 1634 1662 & qair, t_qair, tbnd_qair, & 1635 & rainf, snowf, t_prec, tbnd_prec, prectime,&1663 & rainf, snowf, t_prec, tbnd_prec, & 1636 1664 & swdown, t_swdown, tbnd_swdown, & 1637 1665 & lwdown, t_lwdown, tbnd_lwdown, & … … 1645 1673 REAL(r_std), INTENT(in) :: time_int(2) !! The time interval over which the forcing is needed. 1646 1674 ! 1675 LOGICAL, INTENT(in) :: first_call_readslab 1647 1676 REAL(r_std), INTENT(out) :: tair(:,:), tairmax(:,:), tairmin(:,:) 1648 1677 REAL(r_std), INTENT(out) :: t_tair(:) … … 1657 1686 REAL(r_std), INTENT(out) :: t_prec(:) 1658 1687 REAL(r_std), INTENT(out) :: tbnd_prec(:,:) 1659 REAL(r_std), INTENT(out) :: prectime(:)1660 1688 ! 1661 1689 REAL(r_std), INTENT(out) :: swdown(:,:) … … 1696 1724 CHARACTER(LEN=80) :: cellmethod 1697 1725 ! 1698 LOGICAL, SAVE :: first_call_readslab=.TRUE.1699 !1700 1726 ALLOCATE(time_tmp(slab_size,nbtax)) 1701 1727 ALLOCATE(rau(nbpoint_loc,slab_size)) … … 1713 1739 IF ( first_call_readslab ) THEN 1714 1740 ! 1715 preciptime(:) = 01716 1741 ! 1717 1742 ! If the first file is only there to provide the last time step of the previous year, we … … 1722 1747 current_offset = MAX(imin(1)-2,1) 1723 1748 ! 1724 first_call_readslab = .FALSE.1725 write(numout, *) "first_call_readslab in forcing_readslab_root"1726 1749 ! 1727 1750 ELSE 1728 1751 ! 1729 ! Put back the cummulated time of rainfall into the global array1730 !1731 preciptime(position_slab(1):position_slab(2)) = preciptime(position_slab(1):position_slab(2)) + &1732 & prectime(1:slab_size)1733 !1734 1752 ! Compute new offset 1735 1753 ! 1736 1754 current_offset = position_slab(2)-2 1737 write(numout, *) "first_call_readslab in forcing_readslab_root 22"1738 1755 ! 1739 1756 ENDIF … … 1930 1947 t_prec(1:slab_size) = time(position_slab(1):position_slab(2), timeid_precip) 1931 1948 tbnd_prec(1:slab_size,:) = time_bounds(position_slab(1):position_slab(2),timeid_precip,:) 1932 prectime(1:slab_size) = preciptime(position_slab(1):position_slab(2))1933 1949 ! 1934 1950 WRITE(*,*) "For SWdown we are using time axis '",TRIM(time_axename(timeid_swdown)),& … … 3253 3269 !_ ============================================================================================================================== 3254 3270 3255 SUBROUTINE forcing_zoomgrid (zoom_lon, zoom_lat, filename, dumpncdf)3271 SUBROUTINE forcing_zoomgrid (zoom_lon, zoom_lat, filename,model_guess, dumpncdf) 3256 3272 ! 3257 3273 ! Get the area to be zoomed and the sizes of arrays we will need. … … 3262 3278 ! 3263 3279 REAL(r_std), DIMENSION(2), INTENT(in) :: zoom_lon, zoom_lat 3264 CHARACTER(LEN=*), INTENT(in) :: filename 3280 CHARACTER(LEN=*), INTENT(in) :: filename, model_guess 3265 3281 LOGICAL, INTENT(in) :: dumpncdf 3266 3282 ! … … 3273 3289 REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lon_val, lat_val 3274 3290 INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: zoom_index 3291 INTEGER(i_std),DIMENSION(2) :: tmp_lon_ind1, tmp_lat_ind1, tmp_lon_ind2, tmp_lat_ind2 3292 INTEGER(i_std) :: iindex_init, jindex_init, iindex_end, jindex_end, i_orig, j_orig 3275 3293 ! 3276 3294 INTEGER(i_std) :: iret, force_id, iv … … 3321 3339 ! Determine the size in x and y of the zoom 3322 3340 ! 3323 IF ( zoom_forcing ) THEN 3341 lon_tmp(:) = 0.0 3342 lat_tmp(:) = 0.0 3343 IF ( zoom_forcing .AND. model_guess .EQ. 'WRF') THEN 3344 tmp_lon_ind1 = MINLOC(ABS(longlo_tmp(:,:)-zoom_lon(1))) 3345 tmp_lat_ind1 = MINLOC(ABS(lat_glo(:,:)-zoom_lat(1))) 3346 3347 tmp_lon_ind2 = MINLOC(ABS(longlo_tmp(:,:)-zoom_lon(2))) 3348 tmp_lat_ind2 = MINLOC(ABS(lat_glo(:,:)-zoom_lat(2))) 3349 3350 iindex_init = tmp_lon_ind1(1) 3351 jindex_init = tmp_lat_ind1(2) 3352 3353 iindex_end= tmp_lon_ind2(1) 3354 jindex_end= tmp_lat_ind2(2) 3355 3356 lat_tmp(jindex_init:jindex_end) = 1.0 3357 lon_tmp(iindex_init:iindex_end) = 1.0 3358 3359 iim_loc = iindex_end - iindex_init +1 3360 jjm_loc = jindex_end - jindex_init +1 3361 3362 ELSE IF ( zoom_forcing .AND. model_guess .EQ. 'regular') THEN 3324 3363 ! 3325 3364 ! Working with the hypothesis it is a regular global grid and bring it back to the -180 to 180 interval … … 3346 3385 jjm_loc = NINT(SUM(lat_tmp)) 3347 3386 ELSE 3387 iindex_init = 1 3388 jindex_init = 1 3389 3390 iindex_end= size(longlo_tmp(:,1)) 3391 jindex_end= size(lat_glo(1,:)) 3392 3348 3393 iim_loc = iim_glo 3349 3394 jjm_loc = jjm_glo … … 3383 3428 ! Now find the correspondance between the zoomed & re-ordered grid and the global one. 3384 3429 ! 3385 DO i=1,iim_loc 3386 DO j=1,jjm_loc 3387 ! 3388 imin=MINLOC(ABS(longlo_tmp(:,1)-lon_val(i))) 3389 jmin=MINLOC(ABS(lat_glo(1,:)-lat_val(j))) 3390 ! 3391 lon_loc(i,j) = longlo_tmp(imin(1),jmin(1)) 3392 lat_loc(i,j) = lat_glo(imin(1),jmin(1)) 3393 mask_loc(i,j) = mask_glo(imin(1),jmin(1)) 3394 ! 3395 zoom_index(i,j,1) = imin(1) 3396 zoom_index(i,j,2) = jmin(1) 3397 ! 3430 IF (model_guess .EQ. 'WRF' .AND. zoom_forcing) THEN 3431 ! if WRF grid and zoomed grid, added a few lines here 3432 mask_loc(:,:) = mask_glo(iindex_init:iindex_end, jindex_init:jindex_end) 3433 lat_loc(:,:) = lat_glo(iindex_init:iindex_end, jindex_init:jindex_end) 3434 lon_loc(:,:) = longlo_tmp(iindex_init:iindex_end, jindex_init:jindex_end) 3435 3436 3437 nbpoint_loc = 0 3438 DO i=1,iim_loc 3439 DO j=1,jjm_loc 3440 zoom_index(i,j,1) = iindex_init + i - 1 3441 zoom_index(i,j,2) = jindex_init + j - 1 3442 ! this is consistent to glogrd_getwrf code 3443 if ( mask_loc(i,j) .GT. 0.5 ) THEN 3444 nbpoint_loc = nbpoint_loc +1 3445 ENDIF 3446 ENDDO 3398 3447 ENDDO 3399 ENDDO 3400 ! 3401 nbpoint_loc = SUM(mask_loc) 3448 3449 ELSE 3450 ! if regular grid, or wrf full grids: 3451 ! the orginal code, not modified 3452 DO i=1,iim_loc 3453 DO j=1,jjm_loc 3454 ! 3455 imin=MINLOC(ABS(longlo_tmp(:,1)-lon_val(i))) 3456 jmin=MINLOC(ABS(lat_glo(1,:)-lat_val(j))) 3457 ! 3458 lon_loc(i,j) = longlo_tmp(imin(1),jmin(1)) 3459 lat_loc(i,j) = lat_glo(imin(1),jmin(1)) 3460 mask_loc(i,j) = mask_glo(imin(1),jmin(1)) 3461 ! 3462 zoom_index(i,j,1) = imin(1) 3463 zoom_index(i,j,2) = jmin(1) 3464 ! 3465 ENDDO 3466 ENDDO 3467 ! 3468 nbpoint_loc = SUM(mask_loc) 3469 3470 ENDIF 3402 3471 3403 3472 … … 3482 3551 ! Treat first the longitudes 3483 3552 ! 3484 DO j=1,jjm_loc 3485 dx = zero 3486 DO i=1,iim_loc-1 3487 dx = dx+ABS(lon_loc(i,j)-lon_loc(i+1,j)) 3553 IF ( model_guess .EQ. "WRF" ) THEN 3554 ! if grid is wrf 3555 ! 3556 DO i=1,iim_loc 3557 DO j=1,jjm_loc 3558 i_orig = iindex_init + i - 1 3559 j_orig = jindex_init + j - 1 3560 3561 ! Corners 3562 CALL grid_tolola(i_orig+0.5, j_orig+0.5, corners_loc(i,j,1,1), corners_loc(i,j,1,2)) 3563 CALL grid_tolola(i_orig+0.5, j_orig-0.5, corners_loc(i,j,2,1), corners_loc(i,j,2,2)) 3564 CALL grid_tolola(i_orig-0.5, j_orig-0.5, corners_loc(i,j,3,1), corners_loc(i,j,3,2)) 3565 CALL grid_tolola(i_orig-0.5, j_orig-0.5, corners_loc(i,j,4,1), corners_loc(i,j,4,2)) 3566 3567 coslat = MAX( COS(lat_glo(i_orig,j_orig) * pi/180. ), mincos ) 3568 dx = ABS(corners_loc(i,j,2,1) - corners_loc(i,j,1,1)) * pi/180. * R_Earth * coslat 3569 dy = ABS(corners_loc(i,j,1,2) - corners_loc(i,j,3,2)) * pi/180. * R_Earth 3570 area_loc(i,j) = dx*dy 3571 ENDDO 3488 3572 ENDDO 3489 dx = dx/(iim_loc-1) 3573 3574 ELSE 3575 ! if regular grid: original code, not modified 3576 ! Treat first the longitudes 3577 ! 3578 DO j=1,jjm_loc 3579 dx = zero 3580 DO i=1,iim_loc-1 3581 dx = dx+ABS(lon_loc(i,j)-lon_loc(i+1,j)) 3582 ENDDO 3583 dx = dx/(iim_loc-1) 3584 DO i=1,iim_loc 3585 corners_loc(i,j,1,1) = lon_loc(i,j)-dx/2.0 3586 corners_loc(i,j,2,1) = lon_loc(i,j)+dx/2.0 3587 corners_loc(i,j,3,1) = lon_loc(i,j)+dx/2.0 3588 corners_loc(i,j,4,1) = lon_loc(i,j)-dx/2.0 3589 ENDDO 3590 ENDDO 3591 ! 3592 ! Now treat the latitudes 3593 ! 3490 3594 DO i=1,iim_loc 3491 corners_loc(i,j,1,1) = lon_loc(i,j)-dx/2.0 3492 corners_loc(i,j,2,1) = lon_loc(i,j)+dx/2.0 3493 corners_loc(i,j,3,1) = lon_loc(i,j)+dx/2.0 3494 corners_loc(i,j,4,1) = lon_loc(i,j)-dx/2.0 3595 dy = zero 3596 DO j=1,jjm_loc-1 3597 dy = dy + ABS(lat_loc(i,j)-lat_loc(i,j+1)) 3598 ENDDO 3599 dy = dy/(jjm_loc-1) 3600 DO j=1,jjm_loc 3601 corners_loc(i,j,1,2) = lat_loc(i,j)+dy/2.0 3602 corners_loc(i,j,2,2) = lat_loc(i,j)+dy/2.0 3603 corners_loc(i,j,3,2) = lat_loc(i,j)-dy/2.0 3604 corners_loc(i,j,4,2) = lat_loc(i,j)-dy/2.0 3605 ENDDO 3495 3606 ENDDO 3496 ENDDO 3497 ! 3498 ! Now treat the latitudes 3499 ! 3500 DO i=1,iim_loc 3501 dy = zero 3502 DO j=1,jjm_loc-1 3503 dy = dy + ABS(lat_loc(i,j)-lat_loc(i,j+1)) 3607 3608 ! 3609 ! Compute the areas of the grid (using the simplification that the grid is regular in lon/lat). 3610 ! 3611 DO i=1,iim_loc 3612 DO j=1,jjm_loc 3613 coslat = MAX( COS(lat_loc(i,j) * pi/180. ), mincos ) 3614 dx = ABS(corners_loc(i,j,2,1) - corners_loc(i,j,1,1)) * pi/180. * R_Earth * coslat 3615 dy = ABS(corners_loc(i,j,1,2) - corners_loc(i,j,3,2)) * pi/180. * R_Earth 3616 area_loc(i,j) = dx*dy 3617 ENDDO 3504 3618 ENDDO 3505 dy = dy/(jjm_loc-1) 3506 DO j=1,jjm_loc 3507 corners_loc(i,j,1,2) = lat_loc(i,j)+dy/2.0 3508 corners_loc(i,j,2,2) = lat_loc(i,j)+dy/2.0 3509 corners_loc(i,j,3,2) = lat_loc(i,j)-dy/2.0 3510 corners_loc(i,j,4,2) = lat_loc(i,j)-dy/2.0 3511 ENDDO 3512 ENDDO 3513 ! 3514 ! Compute the areas of the grid (using the simplification that the grid is regular in lon/lat). 3515 ! 3516 DO i=1,iim_loc 3517 DO j=1,jjm_loc 3518 coslat = MAX( COS(lat_loc(i,j) * pi/180. ), mincos ) 3519 dx = ABS(corners_loc(i,j,2,1) - corners_loc(i,j,1,1)) * pi/180. * R_Earth * coslat 3520 dy = ABS(corners_loc(i,j,1,2) - corners_loc(i,j,3,2)) * pi/180. * R_Earth 3521 area_loc(i,j) = dx*dy 3522 ENDDO 3523 ENDDO 3619 3620 ENDIF 3621 3524 3622 ! 3525 3623 ! If requested we read a variable, zoomin and dump the result into a test file. … … 4011 4109 ALLOCATE(time(nb_forcing_steps, nbtax*nbtmethods), time_bounds(nb_forcing_steps,nbtax*nbtmethods,2)) 4012 4110 ALLOCATE(time_axename(nbtax*nbtmethods), time_cellmethod(nbtax*nbtmethods)) 4013 ALLOCATE(preciptime(nb_forcing_steps))4014 4111 ALLOCATE(time_sourcefile(nb_forcing_steps)) 4015 4112 ALLOCATE(time_id(nb_forcing_steps, nbtax)) … … 4262 4359 ENDDO 4263 4360 ! 4264 ! Set to zero the variable for the cummulated time for rainfall4265 !4266 preciptime(:) = zero4267 4361 ! 4268 4362 ! Keep the very first date of the time axis for future use
Note: See TracChangeset
for help on using the changeset viewer.
