source: trunk/LMDZ.MARS/libf/phymars/surfini.F @ 1242

Last change on this file since 1242 was 1226, checked in by aslmd, 11 years ago

LMDZ.MARS : Replaced comcstfi and planete includes by modules.
File size: 20.1 KB
Line 
1      SUBROUTINE surfini(ngrid,piceco2,qsurf,psolaralb)
2   ! to use  'getin'
3      USE ioipsl_getincom
4      use netcdf
5      use tracer_mod, only: nqmx, noms
6      use comgeomfi_h, only: long, lati
7      use surfdat_h, only: watercaptag, frost_albedo_threshold,
8     &                     albedo_h2o_ice, inert_h2o_ice, albedodat,
9     &                     albedice, dryness
10#ifndef MESOSCALE
11      use mod_grid_phy_lmdz, only : klon_glo ! # of physics point on full grid
12      use mod_phys_lmdz_para, only : is_master, gather, scatter
13#endif
14      USE comcstfi_h
15      IMPLICIT NONE
16c=======================================================================
17c
18c   creation des calottes pour l'etat initial
19c
20c=======================================================================
21c-----------------------------------------------------------------------
22c   Declarations:
23c   -------------
24#include "dimensions.h"
25!#include "dimphys.h"
26!#include "surfdat.h"
27#include "callkeys.h"
28!#include "tracer.h"
29!#include "comgeomfi.h"
30
31#include "datafile.h"
32
33      integer,intent(in) :: ngrid ! number of atmospheric columns
34      real,intent(in) :: piceco2(ngrid) ! CO2 ice thickness
35      real,intent(inout) :: qsurf(ngrid,nqmx) ! tracer on surface (kg/m2)
36      real,intent(out) :: psolaralb(ngrid,2) ! albedo
37
38      INTEGER ig,icap,iq,alternate
39      REAL icedryness ! ice dryness
40     
41      ! longwatercaptag is watercaptag. Trick for some compilers
42      LOGICAL, DIMENSION(100000) :: longwatercaptag
43     
44! There are 3 different modes for ice distribution:
45! icelocationmode = 1 ---> based on data from surface.nc
46! icelocationmode = 2 ---> directly predefined for GCM resolutions 32x24 or 64x48
47! icelocationmode = 3 ---> based on logical relations for latitude and longitude
48! For visualisation : > /u/tnalmd/bin/watercaps gcm_txt_output_file
49      INTEGER,SAVE :: icelocationmode = 2
50       
51       
52      !in case icelocationmode == 1
53      INTEGER i,j
54      INTEGER     imd,jmd
55      PARAMETER   (imd=360,jmd=180)
56      REAL        zdata(imd,jmd)
57      REAL        zelat,zelon
58
59#ifndef MESOSCALE
60      INTEGER nb_ice(klon_glo,2)   ! number of counts | detected ice for GCM grid
61#endif
62      INTEGER latice(jjm,2),lonice (iim,2) ! number of counts | detected ice along lat & lon axis
63
64      REAL step,count,ratiolat
65
66      INTEGER   ierr,nid,nvarid
67     
68      REAL,SAVE :: min_icevalue = 500.
69      character(len=50) :: string = 'thermal'
70     
71      character (len=100) :: zedatafile
72
73#ifdef MESOSCALE
74
75      do ig=1,ngrid
76
77         !write(*,*) "all qsurf to zero. dirty."
78         do iq=1,nqmx
79         qsurf(ig,iq)=0.  !! on jette les inputs GCM
80                          !! on regle juste watercaptag
81                          !! il faudrait garder les inputs GCM
82                          !! si elles sont consequentes
83         enddo
84         if ( ( lati(ig)*180./pi .gt. 70. ) .and.
85     .        ( albedodat(ig) .ge. 0.26   ) )  then
86                 write(*,*)"outlier ",ig
87                 watercaptag(ig)  = .true.
88                 dryness(ig)      = 1.
89                 albedodat(ig)    = albedo_h2o_ice  !! pour output
90         else
91                 watercaptag(ig)  = .false.
92                 dryness(ig)      = 1.
93         endif
94
95      enddo
96#endif
97! problem with nested precompiling flags
98
99#ifndef MESOSCALE
100      ! to handle parallel cases
101#if CPP_PARA
102      logical watercaptag_glo(klon_glo)
103      real dryness_glo(klon_glo)
104      real lati_glo(klon_glo)
105      real long_glo(klon_glo)
106#else
107      logical watercaptag_glo(ngrid)
108      real dryness_glo(ngrid)
109      real lati_glo(ngrid)
110      real long_glo(ngrid)
111#endif
112#endif
113
114#ifndef MESOSCALE
115
116c
117c=======================================================================
118! Initialize watercaptag (default is false)
119      watercaptag_glo(:)=.false.
120
121c     water ice outliers
122c     ------------------------------------------
123
124      IF ((water) .and. (caps)) THEN
125     
126c Perennial H20 north cap defined by watercaptag=true (allows surface to be
127c hollowed by sublimation in vdifc).
128
129c We might not want albedodat to be modified because it is used to write
130c restart files. Instead, albedo is directly modified when needed (i.e.
131c if we have watercaptag and no co2 ice), below and in albedocaps.F90
132
133c       "Dryness coefficient" controlling the evaporation and
134c        sublimation from the ground water ice (close to 1)
135c        HERE, the goal is to correct for the fact
136c        that the simulated permanent water ice polar caps
137c        is larger than the actual cap and the atmospheric
138c        opacity not always realistic.
139
140         alternate = 0
141         
142         if (ngrid .ne. 1) then
143           watercaptag(:) = .false.
144           longwatercaptag(:) = .false.
145         endif
146         
147         write(*,*) "surfini: Ice dryness ?"
148         icedryness=1. ! default value
149         call getin("icedryness",icedryness)
150         write(*,*) "surfini: icedryness = ",icedryness
151         dryness (:) = icedryness
152         
153      ! To be able to run in parallel, we work on the full grid
154      ! and dispatch results afterwards
155
156      ! start by geting latitudes and logitudes on full grid
157      ! (in serial mode, this is just a copy)
158      call gather(lati,lati_glo)
159      call gather(long,long_glo)
160
161      if (is_master) then
162
163        IF (ngrid .eq. 1) THEN ! special case for 1d --> do nothing
164     
165         print*, 'ngrid = 1, do no put ice caps in surfini.F'
166
167        ELSE IF (icelocationmode .eq. 1) THEN
168     
169         print*,'Surfini: ice caps defined from surface.nc'
170           
171! This method detects ice as gridded value above min_icevalue in the field "string" from surface.nc
172! Typically, it is for thermal inertia above 500 tiu.
173! Two conditions are verified:
174! 1. GCM ice caps are defined such as area is conserved for a given latitude
175! (the approximation is that all points within the GCM latitude resolution have the same area).
176! 2. caps are placed to fill the GCM points with the most detected ice first.
177     
178
179           
180         zedatafile = trim(datafile)
181 
182       
183         ierr=nf90_open(trim(zedatafile)//'/surface.nc',
184     &   NF90_NOWRITE,nid)
185     
186         IF (ierr.NE.nf90_noerr) THEN
187       write(*,*)'Error : cannot open file surface.nc '
188       write(*,*)'(in phymars/surfini.F)'
189       write(*,*)'It should be in :',trim(zedatafile),'/'
190       write(*,*)'1) You can set this path in the callphys.def file:'
191       write(*,*)'   datadir=/path/to/the/datafiles'
192       write(*,*)'2) If necessary, surface.nc (and other datafiles)'
193       write(*,*)'   can be obtained online on:'
194       write(*,*)' http://www.lmd.jussieu.fr/~forget/datagcm/datafile'
195       CALL ABORT
196         ENDIF
197     
198     
199         ierr=nf90_inq_varid(nid, string, nvarid)
200         if (ierr.ne.nf90_noerr) then
201          write(*,*) 'surfini error, cannot find ',trim(string)
202          write(*,*) ' in file ',trim(zedatafile),'/surface.nc'
203          write(*,*)trim(nf90_strerror(ierr))
204          stop
205         endif
206
207         ierr=nf90_get_var(nid, nvarid, zdata)
208
209         if (ierr.ne.nf90_noerr) then
210          write(*,*) 'surfini: error failed loading ',trim(string)
211          write(*,*)trim(nf90_strerror(ierr))
212          stop
213         endif
214 
215                     
216         ierr=nf90_close(nid)
217 
218
219         nb_ice(:,1) = 1 ! default: there is no ice
220         latice(:,1) = 1
221         lonice(:,1) = 1
222         nb_ice(:,2) = 0
223         latice(:,2) = 0
224         lonice(:,2) = 0
225         !print*,'jjm,iim',jjm,iim ! jjm =  nb lati , iim = nb longi
226
227         ! loop over the GCM grid - except for poles (ig=1 and ngrid)
228         do ig=2,klon_glo-1
229     
230          ! loop over the surface file grid     
231          do i=1,imd
232           do j=1,jmd
233             zelon = i - 180.
234             zelat = 90. - j
235            if ((abs(lati_glo(ig)*180./pi-zelat).le.90./real(jjm)) .and.
236     &        (abs(long_glo(ig)*180./pi-zelon).le.180./real(iim))) then
237              ! count all points in that GCM grid point
238              nb_ice(ig,1) = nb_ice(ig,1) + 1
239              if (zdata(i,j) > min_icevalue)
240                 ! count all detected points in that GCM grid point
241     &           nb_ice(ig,2) = nb_ice(ig,2) + 1
242             endif
243           enddo
244          enddo 
245
246        ! projection of nb_ice on GCM lat and lon axes
247          latice(1+(ig-2)/iim,:) =
248     &     latice(1+(ig-2)/iim,:) + nb_ice(ig,:)
249          lonice(1+mod(ig-2,iim),:) =
250     &     lonice(1+mod(ig-2,iim),:) + nb_ice(ig,:) ! lonice is USELESS ...
251
252         enddo ! of do ig=2,klon_glo-1
253     
254
255     
256         ! special case for poles
257         nb_ice(1,2)   = 1  ! ice prescribed on north pole
258         latice(1,:)   = nb_ice(1,:)
259         lonice(1,:)   = nb_ice(1,:)
260         latice(jjm,:) = nb_ice(ngrid,:)
261         lonice(iim,:) = nb_ice(ngrid,:)
262     
263     
264!      print*, 'latice TOT', latice(:,1)
265!      print*, 'latice FOUND', latice(:,2)
266!      print*, 'lonice TOT', lonice(:,1)
267!      print*, 'lonice FOUND', lonice(:,2)
268     
269!      print*, 'lat ratio', int(real(latice(:,2))/real(latice(:,1))*iim)
270!      print*, 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))*jjm)
271     
272!      print*,''
273!      print*,'sum lat', sum(latice(:,1)), sum(lonice(:,1))
274!      print*,'sum lon', sum(latice(:,2)), sum(lonice(:,2))
275     
276   
277         ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
278         do i=1,jjm/2
279          step  = 1. ! threshold to add ice cap
280          count = 0. ! number of ice GCM caps at this latitude
281          ! ratiolat is the ratio of area covered by ice within this GCM latitude range
282          ratiolat  = real(latice(i,2))/real(latice(i,1))
283          !print*,'i',i,(i-1)*iim+2,i*iim+1
284     
285          ! put ice caps while there is not enough ice,
286          ! as long as the threshold is above 20%
287          do while ( (count .le. ratiolat*iim ) .and. (step .ge. 0.2))
288           count = 0.
289           ! loop over GCM longitudes
290           do j=1,iim
291            ! if the detected ice ratio in the GCM grid point
292            ! is more than 'step', then add ice
293            if (real(nb_ice((i-1)*iim+1+j,2))
294     &        / real(nb_ice((i-1)*iim+1+j,1)) .ge. step) then
295                  watercaptag_glo((i-1)*iim+1+j) = .true.
296                  count = count + 1
297            endif
298           enddo ! of do j=1,iim
299           !print*, 'step',step,count,ratiolat*iim
300           step = step - 0.01
301          enddo ! of do while
302          !print*, 'step',step,count,ratiolat*iim
303
304         enddo ! of do i=1,jjm/2
305           
306
307        ELSE IF (icelocationmode .eq. 2) THEN
308     
309         print*,'Surfini: predefined ice caps'
310     
311         if ((iim .eq. 32) .and. (jjm .eq. 24)) then ! 32x24
312           
313          print*,'water ice caps distribution for 32x24 resolution'
314          longwatercaptag(1:9)    = .true. ! central cap - core
315          longwatercaptag(26:33)  = .true. ! central cap
316          longwatercaptag(1:33)  = .true. ! central cap
317          longwatercaptag(56)  = .true. ! central cap
318          longwatercaptag(58)  = .true. ! central cap
319          longwatercaptag(60)  = .true. ! central cap
320          longwatercaptag(62)  = .true. ! central cap
321          longwatercaptag(64)  = .true. ! central cap
322!---------------------   OUTLIERS  ----------------------------
323
324         else if ((iim .eq. 64) .and. (jjm .eq. 48)) then ! 64x48
325
326          print*,'water ice caps distribution for 64x48 resolution'
327          longwatercaptag(1:65)   = .true. ! central cap - core
328          longwatercaptag(75:85)  = .true. ! central cap
329          longwatercaptag(93:114) = .true. ! central cap
330!---------------------   OUTLIERS  ----------------------------
331          if (.true.) then
332          longwatercaptag(136)    = .true. ! outlier, lat = 78.75
333          longwatercaptag(138)    = .true. ! outlier, lat = 78.75
334          longwatercaptag(140)    = .true. ! outlier, lat = 78.75
335          longwatercaptag(142)    = .true. ! outlier, lat = 78.75
336          longwatercaptag(161)    = .true. ! outlier, lat = 78.75
337          longwatercaptag(163)    = .true. ! outlier, lat = 78.75
338          longwatercaptag(165)    = .true. ! outlier, lat = 78.75
339          longwatercaptag(183)    = .true. ! outlier, lat = 78.75
340          longwatercaptag(185)    = .true. ! outlier, lat = 78.75
341          longwatercaptag(187)    = .true. ! outlier, lat = 78.75
342          longwatercaptag(189)    = .true. ! outlier, lat = 78.75
343          longwatercaptag(191)    = .true. ! outlier, lat = 78.75
344          longwatercaptag(193)    = .true. ! outlier, lat = 78.75
345          longwatercaptag(194)    = .true. ! outlier, lat = 75
346          longwatercaptag(203)    = .true. ! outlier, lat = 75
347          longwatercaptag(207)    = .true. ! outlier, lat = 75
348          longwatercaptag(244)    = .true. ! outlier, lat = 75
349          longwatercaptag(246)    = .true. ! outlier, lat = 75
350          longwatercaptag(250)    = .true. ! outlier, lat = 75
351          longwatercaptag(252)    = .true. ! outlier, lat = 75
352          longwatercaptag(254)    = .true. ! outlier, lat = 75
353          longwatercaptag(256)    = .true. ! outlier, lat = 75
354          endif
355!--------------------------------------------------------------       
356
357
358           
359         else if (klon_glo .ne. 1) then
360       
361          print*,'No predefined ice location for this resolution :',
362     &           iim,jjm
363          print*,'Please change icelocationmode in surfini.F'
364          print*,'Or add some new definitions ...'
365          call abort
366         
367         endif
368
369         do ig=1,klon_glo
370          if (longwatercaptag(ig)) watercaptag_glo(ig) = .true.
371         enddo
372
373
374        ELSE IF (icelocationmode .eq. 3) THEN
375     
376         print*,'Surfini: ice caps defined by lat and lon values'
377
378         do ig=1,klon_glo
379         
380c-------- Towards olympia planitia water caps -----------
381c--------------------------------------------------------
382
383          if ( ( ( lati_glo(ig)*180./pi .ge. 77.  ) .and. ! cap #2
384     .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
385     .           ( long_glo(ig)*180./pi .ge. 110. ) .and.
386     .           ( long_glo(ig)*180./pi .le. 181. ) )
387     .         .or.
388
389     .         ( ( lati_glo(ig)*180./pi .ge. 75.  ) .and. ! cap #4 (Korolev crater)
390     .           ( lati_glo(ig)*180./pi .le. 76.  ) .and.
391     .           ( long_glo(ig)*180./pi .ge. 150. ) .and.
392     .           ( long_glo(ig)*180./pi .le. 168. ) )
393     .         .or.
394     .         ( ( lati_glo(ig)*180./pi .ge. 77 ) .and. ! cap #5
395     .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
396     .           ( long_glo(ig)*180./pi .ge. -150.) .and.
397     .           ( long_glo(ig)*180./pi .le. -110.) ) )
398     .         then
399             
400               if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
401              !    watercaptag(ig)=.true.
402                  alternate = 1
403               else
404                  alternate = 0
405               endif !end if alternate = 0
406               
407          endif
408
409c----------- Opposite olympia planitia water cap --------
410c--------------------------------------------------------
411
412          if ( ( ( lati_glo(ig)*180./pi     .ge.  80 ) .and.
413     .         ( lati_glo(ig)*180./pi     .le.  84 ) )
414     .         .and.
415     .       ( ( long_glo(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
416     .         ( long_glo(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
417!     .     ( ( ( long_glo(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
418!     .         ( long_glo(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
419!     .       ( ( long_glo(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
420!     .         ( long_glo(ig)*180./pi .le. -70. ) ) ) ) then  !!! 64x48
421        !   watercaptag_glo(ig)=.true.
422          endif
423
424
425c -------------------- Central cap ----------------------
426c--------------------------------------------------------
427
428          if (abs(lati_glo(ig)*180./pi).gt.80)
429     .          watercaptag_glo(ig)=.true.
430           
431c--------------------------------------------------------
432c--------------------------------------------------------
433         end do ! of (klon_glo)
434
435
436        ELSE
437     
438         print*, 'In surfini.F, icelocationmode is ', icelocationmode
439         print*, 'It should be 1, 2 or 3.'
440         call abort
441
442        ENDIF ! of if (icelocation)
443       
444       
445        ! print caps locations - useful for plots too
446        print*,'surfini: latitude | longitude | ig'
447        do ig=1,klon_glo
448          dryness_glo(ig) = icedryness
449
450          if (watercaptag_glo(ig)) then
451             print*,'surfini: ice water cap', lati_glo(ig)*180./pi,
452     &              long_glo(ig)*180./pi, ig
453          endif
454        enddo
455       
456       endif !of if (is_master)
457       
458       ! Now scatter fields watercaptag and dryness from master to all
459       ! (is just a plain copy in serial mode)
460       call scatter(dryness_glo,dryness)
461       call scatter(watercaptag_glo,watercaptag)
462       
463! end of #else of #ifdef MESOSCALE
464       ENDIF ! (caps & water)
465#endif       
466
467c ===============================================================
468c      INITIAL ALBEDO
469c ===============================================================
470
471         write(*,*)"surfini: water frost thickness",
472     s     frost_albedo_threshold
473         write(*,*)"surfini: water ice albedo:", albedo_h2o_ice
474         write(*,*)"surfini: water ice TI:", inert_h2o_ice
475
476c        To start with : Initial albedo = observed dataset
477c        -------------------------------------------------
478         DO ig=1,ngrid
479              psolaralb(ig,1)=albedodat(ig)
480              psolaralb(ig,2)=albedodat(ig)
481         END DO
482         PRINT*,'surfini: minimum albedo without water caps',
483     &          minval(albedodat)
484         PRINT*,'surfini: maximum albedo without water caps',
485     &          maxval(albedodat)
486
487c        initial albedo if permanent H2O ice is present
488c        ------------------------------------------------
489         IF ((water) .and. (caps)) THEN
490           DO ig=1,ngrid
491            IF (watercaptag(ig)) THEN
492              psolaralb(ig,1) = albedo_h2o_ice
493              psolaralb(ig,2) = albedo_h2o_ice
494            ENDIF
495           END DO
496           PRINT*,'surfini: minimum albedo with water caps',
497     &            minval(albedodat)
498           PRINT*,'surfini: maximum albedo with water caps',
499     &            maxval(albedodat)
500
501         ENDIF
502
503c      changing initial albedo if CO2 ice is present
504c      -------------------------------------------
505
506       DO ig=1,ngrid
507         IF (piceco2(ig) .GT. 0.) THEN
508             IF(lati(ig).LT. 0.) THEN
509                icap=2 ! Southern hemisphere
510             ELSE
511                icap=1 ! Northern hemisphere
512             ENDIF
513             psolaralb(ig,1) = albedice(icap)
514             psolaralb(ig,2) = albedice(icap)   
515         END IF
516       END DO
517
518c      changing initial albedo if water ice frost is present
519c      -------------------------------------------
520       IF (water) THEN
521          do iq=1,nqmx
522c          if there is frost and surface albedo is set to albedo_h2o_ice
523           if(noms(iq).eq."h2o_ice") then
524             do ig=1,ngrid
525             
526              if ((watercaptag(ig).eqv..false.)
527     &     .and. (qsurf(ig,iq).lt.-frost_albedo_threshold)) then
528              print*, ''
529              print*, '!!! PROBLEM in SURFINI !!!!'
530              print*, 'FOUND NEGATIVE SURFACE ICE VALUE WHERE
531     & WATERCAPTAG IS FALSE'
532              print*, ''
533              print*, 'ig,qsurf,threshold' ,
534     &         ig, qsurf(ig,iq), -frost_albedo_threshold
535              print*, ''
536              print*, '1) Check h2o_ice in startfi and ice
537     & distribution in surfini'
538              print*, '2) Use ini_h2osurf option in newstart'
539              print*, ''
540#ifndef MESOINI
541              CALL ABORT
542#else
543              qsurf(ig,iq) = 0.
544#endif
545            endif
546             
547              if ((piceco2(ig) .eq. 0.).and.
548     &          (qsurf(ig,iq).gt.frost_albedo_threshold)) then
549                     psolaralb(ig,1) = albedo_h2o_ice
550                     psolaralb(ig,2) = albedo_h2o_ice
551c                     PRINT*,'surfini.F frost',
552c     &                  lati(ig)*180./pi, long(ig)*180./pi
553               endif
554              enddo
555           endif
556          end do
557          PRINT*,'surfini: minimum albedo with frost and co2',
558     &            minval(albedodat)
559          PRINT*,'surfini: maximum albedo with frost and co2',
560     &            maxval(albedodat)
561       ELSE
562         watercaptag(:) = .false.
563       END IF ! OF IF(water)
564         
565      RETURN
566      END
Note: See TracBrowser for help on using the repository browser.